Merge branch 'doubled-cell' into rc175: a doubled-cell hypothesis for a half-translation pseudo-symmetry that breaks down with resolution

Narrow, guarded: the doubled pass runs only where the added class grows with resolution and is resolved; the doubled
cell is adopted only where one pseudo-translation explains the added class in every shell and the point group does not
rise. Fixes the one remaining lattice failure (half cell -> deposited cell, R-free 0.574 -> 0.275); p.hkl identical on
the ci tier, the supercell gate set, 78 further supercell-candidate sets and all private sets.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3

docs/RUGNUX_REPORT.md

@@ -595,30 +595,47 @@ On rotation data, section 4 reports whether the lattice the run adopted has inte
 positions it does not index. On 60 frames spread over the sweep, after each frame's own integration, the
 lattice is predicted doubled along all three primitive axes and integrated to 3 Å; the reflections split
 into eight parity classes of h, k and l, of which `0 0 0` is the lattice itself and each of the other
-seven is one index-2 superstructure. Nothing is decided on it: whether a superstructure belongs in the
-cell is as much the depositor's call as the data's, and several crystals whose accepted cell is the
-sub-cell carry one.
+seven is one index-2 superstructure. The class read there is also followed out to the run's own limit,
+counting only half-integer reflections that no lattice reflection on the same frame lies within the
+background radius of, each resolution shell against the median of the six other classes - which a single
+doubling leaves empty, so they measure what sits at a half-integer node without being a reflection. The
+log prints that table (`Supercell probe to the limit`).
 
 **`SUPERCELL_CLASS=`** is the parity class with the most intensity over 20–3 Å.
 **`SUPERCELL_OCCUPANCY_PCT=`** is its mean intensity against the lattice's own reflections on the same
 frames, and **`SUPERCELL_ROCK_PCT=`** (± **`SUPERCELL_ROCK_SE_PCT=`**) the part of it that follows the
 partiality the way a Bragg reflection does, from a fit I = a + b p over the class, on the same scale.
 **`SUPERCELL_I_OVER_SIGMA=`** is its mean I/σ. A class near zero on both is empty. One that is occupied
-and rocks is a superstructure whose reflections this run did not integrate; **`SUPERCELL_DOUBLED_CELL=`**
-is the Niggli-reduced cell the lattice would double to, to give with `-C` to process on it. One that is
-occupied but hardly rocks is diffuse or disordered intensity rather than Bragg reflections.
+and rocks holds reflections of a superstructure this cell does not index; **`SUPERCELL_DOUBLED_CELL=`**
+is the Niggli-reduced cell the lattice doubles to. One that is occupied but hardly rocks is diffuse or
+disordered intensity rather than Bragg reflections.
 
 **`SUPERCELL_POSSIBLE=`** is `TRUE` where the class is measured (`SUPERCELL_I_OVER_SIGMA=` at least
 0.5) and part of it rocks like Bragg reflections (`SUPERCELL_ROCK_PCT=` at least 2 %, three standard
-errors clear of zero). It is advice to check, not a finding: the same numbers come from a real doubled
-cell and from a correct cell with weak ordered intensity between its reflections - or with further
-lattice domains whose spots land on the half-integer positions - and which of the two a structure is
-decided by refinement. Process both settings - the run's cell, and `SUPERCELL_DOUBLED_CELL=` given with
-`-C` - and compare them there. A `TRUE` is also a warning under the `SUPERCELL_POSSIBLE` flag, worded as a
-prompt to check: on a battery of rotation data sets most of the crystals it named refine normally in the
-sub-cell (rocking parts of 2-16 %), but a real doubled cell can read inside that range too. The summary's
-`Supercell` line carries the same advice, and says where further lattice domains were found that may put
-spots on the half-integer positions.
+errors clear of zero). It says the doubled cell is there in the data, not that it is the cell to refine
+in: a crystal whose two halves are nearly alike (a pseudo-translation) is described by either cell, and
+spots of further lattice domains can land on the half-integer positions too. A `TRUE` is also a warning
+under the `SUPERCELL_POSSIBLE` flag, and the summary's `Supercell` line carries the same advice.
+
+**Both cells merged.** Where the class stands clear of the six empty ones out to the limit (in at least
+three shells where it is measured above its noise and most of its reflections are resolved from the
+lattice's, at more than twice the strongest of them in standard errors) and
+grows significantly with resolution to a quarter of the lattice's intensity or more in the finest shell
+both are measured in - the half cell holds at low resolution and not at high, the shape of a
+pseudo-translation breaking down - the run integrates and merges the doubled cell as well, on the same
+orientation and geometry, and decides with both merges in hand. The size of the class does not decide
+it - crystals deposited in the half cell carry classes as strong as one deposited in the doubled cell.
+The doubled cell is the answer only where, in the doubled merge, one translation - half a cell plus a
+small offset, fitted over the class - modulates at least half of the added class in every resolution
+shell, so the two halves of the cell are one content displaced - and where the doubled merge is not more
+symmetric than the half cell's, which a superstructure of it cannot be. A class that is there whatever the
+direction (copies that differ in content, twin domains, neighbouring reflections) keeps the half cell.
+**`DOUBLED_CELL_ADOPTED=`** says which was written as the answer and **`DOUBLED_CELL_DECISION=`** why;
+**`DOUBLED_CELL_OTHER_FILE=`** names the other hypothesis's merged file (`<prefix>_doubled.mtz`, or
+`<prefix>_halfcell.mtz` where the doubled cell was adopted, with the rest of its files under the same
+prefix), and **`DOUBLED_CELL_OTHER_CELL=`**, **`DOUBLED_CELL_OTHER_SPACE_GROUP=`**,
+**`DOUBLED_CELL_OTHER_D_MIN=`**, **`DOUBLED_CELL_OTHER_CC_HALF=`** and, with `--model`,
+**`DOUBLED_CELL_OTHER_R_FREE=`** describe it. Compare the two in refinement.
 
 ## Further lattices
 

image_analysis/IndexAndRefine.cpp

@@ -5,6 +5,7 @@
 #include <cmath>
 #include <limits>
 #include <cstdlib>
+#include <unordered_map>
 #include "IndexAndRefine.h"
 
 #include "bragg_integration/CalcISigma.h"
@@ -633,6 +634,15 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
         ProbeSupercellFrame(msg, prediction, probe_integrate, outcome, settings_prediction);
 }
 
+CrystalLattice DoubledLattice(const CrystalLattice &primitive, const std::array<int, 3> &parity) {
+    const Coord v[3] = {primitive.Vec0(), primitive.Vec1(), primitive.Vec2()};
+    const int i = parity[0] ? 0 : parity[1] ? 1 : 2;
+    Coord t[3];
+    for (int j = 0; j < 3; j++)
+        t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast<float>(parity[j]);
+    return {t[0], t[1], t[2]};
+}
+
 SupercellProbeFit FitSupercellProbe(const SupercellProbeClass &c) {
     SupercellProbeFit f;
     if (c.n < 3)
@@ -707,8 +717,69 @@ void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction
         c.sum_pi += r.partiality * r.I;
         c.sum_ii += static_cast<double>(r.I) * r.I;
     }
+
+    // The same supercell again, to the run's own limit (SupercellTrend). Its own prediction, so the band
+    // above integrates exactly what it always has.
+    BraggPredictionSettings s_trend = settings;
+    s_trend.centering = 'P';
+    ApplyPredictionRange(s_trend, experiment, super);
+    const auto ntrend = prediction.Calc(outcome.experiment, super, s_trend);
+    const std::vector<Reflection> &predicted = prediction.GetReflections();
+
+    // The lattice's own reflections on this frame, binned on a grid of r3 squares so that every one
+    // within r3 of a point is in the 3x3 squares around it.
+    const float r3 = experiment.GetBraggIntegrationSettings().GetR3();
+    const auto square = [r3](float x, float y) {
+        return (static_cast<int64_t>(std::floor(x / r3)) << 32)
+               | static_cast<uint32_t>(static_cast<int32_t>(std::floor(y / r3)));
+    };
+    std::unordered_map<int64_t, std::vector<std::pair<float, float>>> lattice_nodes;
+    for (int i = 0; i < ntrend; i++)
+        if (!((predicted[i].h | predicted[i].k | predicted[i].l) & 1))
+            lattice_nodes[square(predicted[i].predicted_x, predicted[i].predicted_y)].emplace_back(
+                    predicted[i].predicted_x, predicted[i].predicted_y);
+    const auto clear_of_lattice = [&](float x, float y) {
+        for (int dx = -1; dx <= 1; dx++)
+            for (int dy = -1; dy <= 1; dy++) {
+                const auto it = lattice_nodes.find(square(x + dx * r3, y + dy * r3));
+                if (it == lattice_nodes.end())
+                    continue;
+                for (const auto &[nx, ny] : it->second)
+                    if ((nx - x) * (nx - x) + (ny - y) * (ny - y) < r3 * r3)
+                        return false;
+            }
+        return true;
+    };
+
+    SupercellTrend trend{};
+    for (const auto &r : integrate(predicted, ntrend, msg.number)) {
+        if (!(r.d <= SUPERCELL_TREND_D_A[0]) || !(r.sigma > 0.0f) || !std::isfinite(r.I))
+            continue;
+        const int c = 4 * (r.h & 1) + 2 * (r.k & 1) + (r.l & 1);
+        if (c != 0 && !clear_of_lattice(r.predicted_x, r.predicted_y))
+            continue;
+        int sh = 0;
+        while (r.d <= SUPERCELL_TREND_D_A[sh + 1])
+            sh++;
+        auto &t = trend[c][sh];
+        t.n++;
+        t.sum_i += r.I;
+        t.sum_i_over_sigma += r.I / r.sigma;
+        t.sum_ii += static_cast<double>(r.I) * r.I;
+    }
+
     const std::unique_lock ul(supercell_probe_mutex_);
-    supercell_probe_by_frame_.insert_or_assign(msg.number, ProbedFrame{frame, prim});
+    supercell_probe_by_frame_.insert_or_assign(msg.number, ProbedFrame{frame, trend, prim});
+}
+
+SupercellTrend IndexAndRefine::GetSupercellTrend() {
+    const std::unique_lock ul(supercell_probe_mutex_);
+    SupercellTrend sum{};
+    for (const auto &[image, f] : supercell_probe_by_frame_)
+        for (int i = 0; i < 8; i++)
+            for (int j = 0; j < SUPERCELL_TREND_SHELLS; j++)
+                sum[i][j] += f.trend[i][j];
+    return sum;
 }
 
 std::optional<IndexAndRefine::IndexingOutcome>

image_analysis/IndexAndRefine.h

@@ -50,6 +50,22 @@ struct SupercellProbeClass {
 // [class][shell], the shells 20-5 A and 5-3 A.
 using SupercellProbe = std::array<std::array<SupercellProbeClass, 2>, 8>;
 
+// The supercell probe's reach to the run's resolution limit: the eight parity classes again, in fixed
+// resolution shells from 20 A down to the prediction limit (SUPERCELL_TREND_D_A gives the shell edges),
+// with a half-integer reflection counted only where no lattice reflection predicted on the same frame
+// lies within the background radius r3 of it - so what is summed is its own intensity, not the tail of
+// a neighbour's. A superstructure from a pseudo-translation is weakest at low resolution, where the two
+// halves of the cell look alike, and grows towards the limit, which the 20-3 A band cannot show.
+inline constexpr float SUPERCELL_TREND_D_A[] = {20.0f, 8.0f, 5.0f, 4.0f, 3.3f, 2.8f, 2.4f, 2.1f, 1.9f, 1.7f,
+                                                1.5f, 1.3f, 1.1f, 0.0f};
+inline constexpr int SUPERCELL_TREND_SHELLS = std::size(SUPERCELL_TREND_D_A) - 1;
+using SupercellTrend = std::array<std::array<SupercellProbeClass, SUPERCELL_TREND_SHELLS>, 8>;
+
+// The index-2 supercell of a primitive lattice whose added reflections are the half-integer class of
+// parity (h,k,l): the vectors t with (h,k,l).t even - twice one axis the class is odd on, and each
+// other axis plus that one where the class is odd on it too. Not reduced.
+CrystalLattice DoubledLattice(const CrystalLattice &primitive, const std::array<int, 3> &parity);
+
 // One parity class of the supercell probe, fitted I = a + b p over partiality p: b (with its standard
 // error) is what rocks like a Bragg reflection, a what does not.
 struct SupercellProbeFit {
@@ -131,6 +147,7 @@ class IndexAndRefine {
     std::vector<int64_t> supercell_probe_frames_;
     struct ProbedFrame {
         SupercellProbe probe;
+        SupercellTrend trend;
         CrystalLattice primitive;
     };
     std::map<int64_t, ProbedFrame> supercell_probe_by_frame_;
@@ -200,6 +217,7 @@ public:
     // The probe's integrations are the caller's to keep out of anything else it counts.
     void ProbeSupercell(std::vector<int64_t> image_numbers);
     [[nodiscard]] SupercellProbe GetSupercellProbe();
+    [[nodiscard]] SupercellTrend GetSupercellTrend();
     // The primitive lattice the parity classes are indexed in (the highest probed frame's; only its metric
     // is meant to be read), so a class can be turned into the cell it would double to.
     [[nodiscard]] std::optional<CrystalLattice> GetSupercellProbePrimitive();

image_analysis/scale_merge/SearchSpaceGroup.cpp

@@ -1676,7 +1676,8 @@ SearchSpaceGroupResult SearchSpaceGroup(
             const double u = std::max(0.0, e) / mean;
             violations += violation;
             sum_u += u;
-            if (!(violation && u >= 1.0))
+            // At its row's mean as MEASURED, not as de-modulated (see the screw zone below).
+            if (!(violation && Esq[i] >= mean))
                 max_u = std::max(max_u, u);
         }
         const int a = static_cast<int>(odd.size());
@@ -2500,7 +2501,15 @@ SearchSpaceGroupResult SearchSpaceGroup(
                 // eight 0k0 are genuinely dead carries one violation at 0.31 of its row - the
                 // mis-measurement the trim exists for, and one that must stay trimmable or a real
                 // 2_1 is lost. Everything below the row mean is trimmed exactly as before.
-                if (!(violation && u >= 1.0))
+                //
+                // "At its row's mean" is read on the MEASURED intensity, not the de-modulated one.
+                // Under a pseudo-translation u is the reflection divided by the depth the modulation
+                // puts on its class, so a reflection at a hundredth of its row reads at its mean once
+                // the class is held at a two-hundredth: measured on a doubled orthorhombic cell, one
+                // 00l at 0.015 of its row (5.9 sigma, beside reflections at 25) became the one member
+                // no trim could touch, and on its own it decided the row against the screw that the
+                // other 26 dead reflections - and the deposited structure factors - carry.
+                if (!(violation && a.e_squared >= mit->second))
                     z.max_u = std::max(z.max_u, u);
                 z.n_control = static_cast<int>(row_present_esq.at(a.row).size());
             }

rugnux/ResultReport.cpp

@@ -1426,10 +1426,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
         }
 
         // ---- index-2 superstructure
-        // Report only, like the harmonic above: on 60 frames the probe integrates the lattice doubled
-        // along all three primitive axes and reads the seven half-integer parity classes. Whether a
-        // class that is there belongs in the cell is the depositor's call as much as the data's -
-        // several crystals whose accepted cell is the sub-cell carry one - so no verdict is printed.
+        // On 60 frames the probe integrates the lattice doubled along all three primitive axes and reads
+        // the seven half-integer parity classes. Whether a class that is there belongs in the cell is the
+        // depositor's call as much as the data's - several crystals whose accepted cell is the sub-cell
+        // carry one - so where the run merged both cells (DOUBLED_CELL_*) the report names both.
         if (result.supercell) {
             const auto &sc = *result.supercell;
             Add(s, KeyText("SUPERCELL_CLASS", fmt::format("{} {} {}", sc.h, sc.k, sc.l)));
@@ -1445,25 +1445,40 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                 Warn(doc, PathologyCode::SUPERCELL_POSSIBLE,
                      fmt::format("Check the cell: half-integer reflections of parity {} {} {} carry Bragg-like intensity "
                                  "({:.0f}% of the lattice's, rocking part {:.1f}%), so the true cell may be twice as "
-                                 "large - or this is weak ordered intensity of a correct cell{}. Process both settings - "
-                                 "this cell and -C \"{:.2f},{:.2f},{:.2f},{:.2f},{:.2f},{:.2f}\" - and compare them in refinement",
+                                 "large - or the crystal's two halves are nearly alike and either cell describes it{}. {}",
                                  sc.h, sc.k, sc.l, sc.occupancy_pct, sc.rock_pct,
-                                 domains ? ", or spots of the further lattice domains found landing on those positions" : "",
-                                 sc.doubled_cell.a, sc.doubled_cell.b, sc.doubled_cell.c,
-                                 sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma));
+                                 domains ? ", or spots of the further lattice domains found land on those positions" : "",
+                                 result.doubled_cell
+                                 ? fmt::format("Both cells were merged; the other one is in {}.mtz - compare them in refinement",
+                                               result.doubled_cell->other_prefix)
+                                 : fmt::format("The doubled cell {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} was not merged",
+                                               sc.doubled_cell.a, sc.doubled_cell.b, sc.doubled_cell.c,
+                                               sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma)));
+            if (result.doubled_cell) {
+                const auto &h = *result.doubled_cell;
+                Add(s, KeyBool("DOUBLED_CELL_ADOPTED", h.adopted));
+                Add(s, KeyText("DOUBLED_CELL_DECISION", h.decision));
+                Add(s, KeyText("DOUBLED_CELL_OTHER_FILE", h.other_prefix + ".mtz"));
+                Add(s, KeyText("DOUBLED_CELL_OTHER_CELL", CellString(h.other_cell)));
+                Add(s, KeyText("DOUBLED_CELL_OTHER_SPACE_GROUP", h.other_space_group));
+                Add(s, KeyReal("DOUBLED_CELL_OTHER_D_MIN", h.other_d_min_A, "{:.2f}"));
+                Add(s, KeyReal("DOUBLED_CELL_OTHER_CC_HALF", h.other_cc_half, "{:.4f}"));
+                if (std::isfinite(h.other_r_free))
+                    Add(s, KeyReal("DOUBLED_CELL_OTHER_R_FREE", h.other_r_free, "{:.4f}"));
+            }
             Add(s, Blank());
             Add(s, Prose(
                     " SUPERCELL_CLASS is the half-integer class (parity of h, k, l in the primitive cell of the\n"
                     " lattice) with the most intensity over 20-3 A, read on 60 frames. OCCUPANCY is its mean\n"
                     " intensity against the lattice's own reflections; ROCK is the part of it that follows the\n"
                     " partiality the way a Bragg reflection does, on the same scale. A class near zero on both is\n"
-                    " empty. One that is occupied and rocks is an index-2 superstructure the lattice does not\n"
-                    " index: its reflections are not integrated. To process on the doubled lattice, give\n"
-                    " SUPERCELL_DOUBLED_CELL with -C. An occupied class that does not rock is diffuse or\n"
-                    " disordered intensity rather than Bragg reflections. SUPERCELL_POSSIBLE is TRUE, and a\n"
-                    " warning to check, where the class is measured and rocks; correct cells with weak ordered\n"
-                    " intensity, or with further lattice domains whose spots land on the half-integer\n"
-                    " positions, read the same, so only refinement of both cells decides."));
+                    " empty. One that is occupied and rocks holds reflections of an index-2 superstructure the\n"
+                    " lattice does not index, and they are not integrated in this cell. An occupied class that\n"
+                    " does not rock is diffuse or disordered intensity rather than Bragg reflections.\n"
+                    " SUPERCELL_POSSIBLE is TRUE, and a warning to check, where the class is measured and rocks.\n"
+                    " It says the doubled cell is there in the data, not that it is the cell to refine in: a\n"
+                    " crystal whose two halves are nearly alike (a pseudo-translation) is described in either,\n"
+                    " and spots of further lattice domains can also land on the half-integer positions."));
             Add(s, Blank());
         }
 
@@ -2783,8 +2798,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                         checked_row("Supercell", {PathologyCode::SUPERCELL_POSSIBLE}, "no indication");
                     else
                         checked_row("Supercell", {PathologyCode::SUPERCELL_POSSIBLE}, fmt::format("POSSIBLE, check - class {} {} {} rocks like Bragg reflections "
-                                                     "(rocking part {:.1f}%); correct cells with weak ordered intensity "
-                                                     "read the same{}; {}", sc.h, sc.k, sc.l, sc.rock_pct,
+                                                     "(rocking part {:.1f}%); a nearly exact pseudo-translation is described "
+                                                     "in either cell{}; {}", sc.h, sc.k, sc.l, sc.rock_pct,
                                                      domains ? ", and the further lattice domains found may put "
                                                                "spots there" : "",
                                                      doubled));

rugnux/Rugnux.cpp

@@ -879,12 +879,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     // at the half cell by every arbiter. The probe reads the question in intensities instead, on the
     // same frames the first pass validates on (see IndexAndRefine::ProbeSupercell).
     //
-    // REPORT ONLY - it decides nothing. Measured on the battery, a weak real class is not yet told apart
-    // from what sits at the half-integer nodes of crystals whose cell is right: a class at a fifth of
-    // the lattice's intensity on a crystal whose axis is truly doubled, against up to a tenth on
-    // correct cells (diffuse scattering and the tails of long-axis neighbours, which do not rock like
-    // Bragg reflections - but neither, on that crystal, does its real class). The line it logs is the
-    // population the decision has to be calibrated on.
+    // It decides nothing in this pass. Read on its own it does not tell a weak real class from what sits
+    // at the half-integer nodes of crystals whose cell is right, so the run carries the doubled cell as a
+    // second hypothesis where the class stands clear out to the limit (DoubledCellHypothesis), and merges
+    // both.
     if (full && indexer && experiment_.IsRotationIndexing()) {
         std::vector<int64_t> frames;
         for (const int o : select_equally_spaced_image_ordinals(images_to_process, 60))

rugnux/Rugnux.h

@@ -458,11 +458,24 @@ struct ProcessResult {
     // The strongest index-2 superstructure class the supercell probe read (20-3 A, 60 frames): which
     // parity class of the primitive cell, its mean intensity and the part that rocks like a Bragg
     // reflection, both against the lattice's own, and the cell the lattice would double to were it
-    // real. Report only - nothing decides on it.
+    // real. What reads it is the doubled-cell hypothesis (Rugnux::DoubledCellHypothesis).
     struct SupercellClass {
         int h = 0, k = 0, l = 0;   // parity of the primitive indices
         double occupancy_pct = 0.0, rock_pct = 0.0, rock_se_pct = 0.0, i_over_sigma = 0.0;
         UnitCell doubled_cell{};
+        // The class out to the run's limit (SupercellTrend), by resolution shell, less the median of the
+        // six other classes: its mean intensity against the lattice's (excess) and in its own standard
+        // errors (z), the largest such z among the six beside it, the lattice's own <I/sigma>, and the
+        // share of the class resolved from the lattice.
+        struct Shell {
+            float d_low_A = 0.0f, d_high_A = 0.0f;
+            int64_t n = 0;
+            double excess = 0.0, z = 0.0, sibling_z = 0.0, lattice_i_over_sigma = 0.0;
+            // The class's reflections clear of the lattice's by r3, per lattice reflection: how much of
+            // the class this geometry resolves at all.
+            double resolved = 0.0;
+        };
+        std::vector<Shell> trend;
     };
     std::optional<SupercellClass> supercell;
     // Rotation: what the spots the crystal's lattice leaves over index to (split-crystal domains, twin
@@ -527,6 +540,23 @@ struct ProcessResult {
     };
     std::optional<TwinSubgroupMerge> twin_subgroup_merge;
 
+    // The doubled-cell hypothesis (Rugnux::DoubledCellHypothesis): where the supercell probe's class
+    // stands clear of the crystal's own null out to the limit, the run integrates and merges in the
+    // doubled cell as well, and one of the two cells is the answer. What the other one is, and where its
+    // merged files are. Unset where the hypothesis was not run.
+    struct DoubledCell {
+        bool adopted = false;              // the answer is the doubled cell, the other one the probed cell
+        std::string decision;              // why, in the words the run reports it in
+        std::string other_prefix;          // <prefix>_halfcell or <prefix>_doubled
+        UnitCell other_cell{};
+        std::string other_space_group;
+        double other_d_min_A = NAN, other_cc_half = NAN, other_r_free = NAN;
+    };
+    std::optional<DoubledCell> doubled_cell;
+    // The merged reflections as written, kept only where the run asks for them (the doubled-cell pass,
+    // whose merge the hypothesis reads). Empty otherwise.
+    std::vector<MergedReflection> merged_reflections;
+
     // What --model found, when one was given. Set even when the validation could not run - its
     // failure_reason then says why - so the report can report the failure rather than leave a run
     // that was asked to validate looking like one that was not. Unset when no model was asked for.
@@ -644,6 +674,36 @@ constexpr float OVERLOAD_LOW_RES_A = 4.0f;
 // reports it in, or an empty string when the refined pass stands.
 std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &pass2);
 
+// The doubled-cell hypothesis (Rugnux::DoubledCellHypothesis), read off the supercell probe's class out
+// to the limit (SupercellClass::trend). Worth merging: the class stands clear of the crystal's own null.
+// Rises: the class grows significantly with resolution, the shape of a pseudo-translation breaking down;
+// empty where it does not, otherwise the reason in the words the run reports it in.
+bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc);
+std::string DoubledCellRises(const ProcessResult::SupercellClass &sc);
+
+// How much of a doubled-cell merge's added class ONE pseudo-translation explains. Two copies related by
+// half a doubled-cell vector plus eps put |F|^2 sin^2(pi h.eps) on the added reflections and cos^2 on
+// the others, so the added class, in units of the others' shell mean, is fitted as A sin^2(pi h.eps) + B:
+// eps over the whole class, A and B shell by shell. B is what is there whatever the direction - copies
+// that differ in content, or intensity that is not the crystal's (twin domains, neighbours) - and
+// `modulated_share` is 1 - B over the shell's mean. `tested` is false where the added class could not
+// be identified in the written cell (no parity class whose sub-lattice is the half cell) or too few of
+// its reflections were measured.
+struct DoubledClassTranslation {
+    struct Shell {
+        float d_low_A = 0.0f, d_high_A = 0.0f;
+        int64_t n = 0;
+        double mean = 0.0, modulated_share = 0.0;
+    };
+    bool tested = false;
+    std::array<int, 3> parity{};
+    std::array<double, 3> eps{};
+    std::vector<Shell> shells;
+};
+DoubledClassTranslation FitDoubledClassTranslation(const std::vector<MergedReflection> &merged,
+                                                   const UnitCell &doubled_cell, const UnitCell &half_cell,
+                                                   char half_centering = 'P');
+
 // Whether two passes hold the same lattice: the same Bravais class and primitive volumes within 2 %.
 bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b);
 
@@ -756,6 +816,10 @@ class Rugnux {
     // Set by Run()'s supercell-collapse guard to force pass-1's full result into the second-pass indexer
     // (RunPipeline calls indexer->ForceRotationIndexerResult when this is set); reset right after that re-run.
     std::optional<RotationIndexerResult> force_rotation_result_;
+    // The indexing result the last full pass integrated on, which the doubled-cell hypothesis doubles.
+    std::optional<RotationIndexerResult> last_rotation_result_;
+    // Keep the written merged reflections on the result (ProcessResult::merged_reflections).
+    bool keep_merged_reflections_ = false;
 
     // The space group the USER fixed (-S on the CLI, the settings panel in the viewer), captured before
     // anything can clear it - a pass that determines a group of its own writes it onto the experiment,
@@ -1051,6 +1115,10 @@ class Rugnux {
     };
     // The refined passes of the rotation two-pass, from the second pass on. Empty where they ran beside
     // the rotation-scale walk and it overturned them.
+    // After the refined passes: where the supercell probe's class stands clear, integrate and merge the
+    // doubled cell too, and return whichever of the two is the answer (ProcessResult::doubled_cell).
+    ProcessResult DoubledCellHypothesis(RugnuxObserver *observer, ProcessResult half,
+                                        const std::string &base_prefix);
     std::optional<ProcessResult> RefinedPasses(RugnuxObserver *observer, ProcessResult pass1,
                                                const RefinedPassStart &start);
     // The rotation-scale walk, run on the run while the refined passes run on a copy of it; `check` is

rugnux/RugnuxImagePass.cpp

@@ -424,9 +424,7 @@ void Rugnux::ImagePass(PipelineLocals &p) {
                         "mean intensity, and the fit I = a + b p against the lattice's own):{}", table);
 
             // The class with the most intensity over 20-3 A goes into the report, with the cell it would
-            // double the lattice to: the vectors t of the primitive lattice with (h,k,l).t even - twice
-            // one axis the class is odd on, and each other axis plus that one where the class is odd
-            // on it too - Niggli-reduced.
+            // double the lattice to (DoubledLattice), Niggli-reduced.
             const auto band = [&](int c) { SupercellProbeClass x = probe[c][0]; x += probe[c][1]; return x; };
             const SupercellProbeFit f0 = FitSupercellProbe(band(0));
             const auto prim = indexer->GetSupercellProbePrimitive();
@@ -436,19 +434,62 @@ void Rugnux::ImagePass(PipelineLocals &p) {
                     best = c;
             if (prim && f0.mean_i > 0.0 && f0.b != 0.0) {
                 const SupercellProbeFit f = FitSupercellProbe(band(best));
-                const int p[3] = {(best >> 2) & 1, (best >> 1) & 1, best & 1};
-                const Coord v[3] = {prim->Vec0(), prim->Vec1(), prim->Vec2()};
-                const int i = p[0] ? 0 : p[1] ? 1 : 2;
-                Coord t[3];
-                for (int j = 0; j < 3; j++)
-                    t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast<float>(p[j]);
+                const std::array<int, 3> p = {(best >> 2) & 1, (best >> 1) & 1, best & 1};
                 ProcessResult::SupercellClass sc;
                 sc.h = p[0]; sc.k = p[1]; sc.l = p[2];
                 sc.occupancy_pct = 100.0 * f.mean_i / f0.mean_i;
                 sc.rock_pct = 100.0 * f.b / f0.b;
                 sc.rock_se_pct = 100.0 * f.b_se / f0.b;
                 sc.i_over_sigma = f.mean_i_over_sigma;
-                sc.doubled_cell = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell();
+                sc.doubled_cell = DoubledLattice(*prim, p).NiggliReduce().GetUnitCell();
+
+                // The same class out to the limit, half-integer reflections clear of the lattice's only
+                // (SupercellTrend). The six other classes, which a single doubling leaves empty, are the
+                // crystal's own null: measured on the same frames through the same integration, they
+                // carry whatever sits at a half-integer node without being a reflection, so the class
+                // is read against their median - in units of the lattice's intensity (excess) and of
+                // its own standard error (z), the strongest sibling's z beside it.
+                const SupercellTrend trend = indexer->GetSupercellTrend();
+                const auto mean = [](const SupercellProbeClass &c) { return c.sum_i / static_cast<double>(c.n); };
+                const auto se = [&](const SupercellProbeClass &c) {
+                    const double m = mean(c);
+                    return std::sqrt(std::max(0.0, c.sum_ii / static_cast<double>(c.n) - m * m)
+                                     / static_cast<double>(c.n));
+                };
+                std::string trend_table;
+                for (int sh = 0; sh < SUPERCELL_TREND_SHELLS; sh++) {
+                    const SupercellProbeClass &l = trend[0][sh], &x = trend[best][sh];
+                    std::vector<int> siblings;
+                    for (int c = 1; c < 8; c++)
+                        if (c != best && trend[c][sh].n > 1)
+                            siblings.push_back(c);
+                    if (l.n == 0 || x.n < 2 || siblings.size() < 6 || l.sum_i <= 0.0)
+                        continue;
+                    std::vector<double> sib_mean;
+                    for (const int c : siblings)
+                        sib_mean.push_back(mean(trend[c][sh]));
+                    std::ranges::sort(sib_mean);
+                    const double null_i = (sib_mean[2] + sib_mean[3]) / 2.0;
+                    ProcessResult::SupercellClass::Shell shell;
+                    shell.d_low_A = SUPERCELL_TREND_D_A[sh];
+                    shell.d_high_A = SUPERCELL_TREND_D_A[sh + 1];
+                    shell.n = x.n;
+                    shell.lattice_i_over_sigma = l.sum_i_over_sigma / static_cast<double>(l.n);
+                    shell.resolved = static_cast<double>(x.n) / static_cast<double>(l.n);
+                    shell.excess = (mean(x) - null_i) / mean(l);
+                    shell.z = (mean(x) - null_i) / se(x);
+                    shell.sibling_z = -INFINITY;
+                    for (const int c : siblings)
+                        shell.sibling_z = std::max(shell.sibling_z, (mean(trend[c][sh]) - null_i) / se(trend[c][sh]));
+                    sc.trend.push_back(shell);
+                    trend_table += fmt::format("\n  {:.1f}-{:.1f} A: lattice n {} <I/s> {:.2f}; class n {} excess {:.3f} "
+                                               "z {:.1f}, other classes z up to {:.1f}", shell.d_low_A, shell.d_high_A,
+                                               l.n, shell.lattice_i_over_sigma, x.n, shell.excess, shell.z,
+                                               shell.sibling_z);
+                }
+                logger.Info("Supercell probe to the limit, class {}{}{} (its reflections clear of the lattice's by r3, "
+                            "less the median of the six other classes: mean intensity against the lattice's own, "
+                            "and in standard errors):{}", p[0], p[1], p[2], trend_table);
                 result.supercell = sc;
             }
         }
@@ -531,6 +572,7 @@ void Rugnux::ImagePass(PipelineLocals &p) {
 
     if (full && !cancelled_) {
         if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) {
+            last_rotation_result_ = *rot;
             end_msg.rotation_lattice = rot->lattice;
             // The refined geometry, not the nominal StartMessage values. Two things about these five
             // numbers that are not obvious from their names. The beam centre is the PONI - the foot of

rugnux/RugnuxPasses.cpp

@@ -19,6 +19,7 @@
 #include <cmath>
 #include <condition_variable>
 #include <deque>
+#include <filesystem>
 #include <fstream>
 #include <cstring>
 #include <functional>
@@ -599,7 +600,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
                 if (!check->overturns) {
                     if (output_error)
                         std::rethrow_exception(output_error);
-                    return *result;
+                    return output->DoubledCellHypothesis(observer, std::move(*result), base_prefix);
                 }
             }
 
@@ -620,7 +621,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
                 logger.Info("Two-pass: the rotation-scale walk changed the spot finding the refined passes "
                             "start from, so they run again after it");
         }
-        return *RefinedPasses(observer, std::move(pass1), start);
+        return DoubledCellHypothesis(observer, *RefinedPasses(observer, std::move(pass1), start), base_prefix);
     }
     return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
 }
@@ -1414,3 +1415,354 @@ std::optional<ProcessResult> Rugnux::RefinedPasses(RugnuxObserver *observer, Pro
     pass2.geometry_not_converged = geometry_not_converged;
     return pass2;
 }
+
+// The supercell probe reads one half-integer class out to the limit (SupercellClass::trend), each
+// shell less the median of the six other classes - which a single doubling leaves empty, so they are the
+// null measured on the same frames, through the same integration and the same neighbour filter.
+//
+// Measured on the forced doubled-cell merges of six crystals whose class the probe reads: a class read
+// this way is a class the doubled merge holds (odd-class half-set CC1/2 0.7-0.98 at 8-4 A on all six),
+// and its SIZE did not tell the crystal deposited in the doubled cell from those deposited in the half
+// cell - the odd class reached 0.4-0.85 of the lattice at high resolution on three of them. Two things
+// about it did, each a significance test rather than a bound:
+//  - its SHAPE with resolution. A pseudo-translation by half a lattice vector leaves the two halves of
+//    the cell alike at low resolution and different at high, so the class grows from nearly nothing
+//    outwards; the half cell is then right at low resolution and wrong at high. Classes that stay a
+//    flat tenth of the lattice at every resolution, or are strongest at low resolution, are not that.
+//  - whether one translation EXPLAINS it (FitDoubledClassTranslation). Copies related by a translation
+//    put the added class where sin^2(pi h.eps) is large and nothing where it is small. On the doubled
+//    deposition it is all modulation (the part left over -0.003 +- 0.005 of the lattice); on the three
+//    rising or strong classes whose depositions kept the half cell 0.08-0.31 is there whatever the
+//    direction, 32-64 standard errors from zero - copies that differ in content, or intensity that is
+//    not the crystal's own.
+namespace {
+    // A shell the class is read in: measured above its noise, against a lattice that is measured too, in a
+    // geometry that resolves most of it from the lattice. Where the added reflections mostly fall within
+    // the background radius of a lattice reflection, what the few clear ones read is not the class, and a
+    // doubled merge integrates the rest on top of their neighbours (measured: a 701 A axis left 26-43 of
+    // ~36000 clear per shell, read them at the lattice's own strength, and its doubled merge passed every
+    // other test).
+    bool DoubledCellShellMeasured(const ProcessResult::SupercellClass::Shell &sh) {
+        return sh.z >= 3.0 && sh.lattice_i_over_sigma >= 1.0 && sh.resolved >= 0.5;
+    }
+}
+
+bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc) {
+    // Clear of the null: in at least three shells the class is measured and stands at more than twice
+    // the strongest of the six empty classes, in standard errors.
+    int clear = 0;
+    for (const auto &sh : sc.trend)
+        clear += DoubledCellShellMeasured(sh) && sh.z > 2.0 * sh.sibling_z;
+    return clear >= 3;
+}
+
+std::string DoubledCellRises(const ProcessResult::SupercellClass &sc) {
+    // From the coarsest shell to the finest one both are measured in, the class has to grow by three
+    // standard errors of the difference, and end at a quarter of the lattice's intensity or more. The
+    // rise alone is not enough: on strong data a class at one or two percent of the lattice rises
+    // significantly (measured: 0.009 to 0.059 at 4 standard errors on a triclinic crystal, and 0.000 to
+    // 0.060 on a cubic one), and a half cell missing a few percent of the intensity still describes the
+    // crystal. The quarter sits far from both sides: the doubled deposition ends at 0.35-0.40, those two
+    // at 0.06.
+    constexpr double DOUBLED_CELL_OUTER_EXCESS = 0.25;
+    if (sc.trend.empty() || !DoubledCellShellMeasured(sc.trend.front()))
+        return {};
+    const auto &low = sc.trend.front();
+    const ProcessResult::SupercellClass::Shell *outer = nullptr;
+    for (const auto &sh : sc.trend)
+        if (DoubledCellShellMeasured(sh))
+            outer = &sh;
+    const double se = std::hypot(low.excess / low.z, outer->excess / outer->z);
+    if (outer == &low || outer->excess - low.excess < 3.0 * se || outer->excess < DOUBLED_CELL_OUTER_EXCESS)
+        return {};
+    return fmt::format("the added class grows from {:.3f} of the lattice's intensity at {:.0f}-{:.0f} A to {:.3f} at "
+                       "{:.1f}-{:.1f} A ({:.0f} standard errors): the half cell holds at low resolution and not at high",
+                       low.excess, low.d_low_A, low.d_high_A, outer->excess, outer->d_low_A, outer->d_high_A,
+                       (outer->excess - low.excess) / se);
+}
+
+DoubledClassTranslation FitDoubledClassTranslation(const std::vector<MergedReflection> &merged,
+                                                   const UnitCell &doubled_cell, const UnitCell &half_cell,
+                                                   char half_centering) {
+    DoubledClassTranslation out;
+    // The added class in the written cell: the parity whose sub-lattice is the half cell. Its real lattice
+    // is the doubled one plus half of (p0 a + p1 b + p2 c), so one vector the class is odd on is replaced
+    // by that half-vector.
+    const CrystalLattice dl(doubled_cell);
+    const Coord v[3] = {dl.Vec0(), dl.Vec1(), dl.Vec2()};
+    const UnitCell want = CrystalLattice(half_cell).ToPrimitive(half_centering).NiggliReduce().GetUnitCell();
+    double best_mismatch = INFINITY;
+    for (int c = 1; c < 8; c++) {
+        const std::array<int, 3> p = {(c >> 2) & 1, (c >> 1) & 1, c & 1};
+        const int i = p[0] ? 0 : p[1] ? 1 : 2;
+        Coord t[3] = {v[0], v[1], v[2]};
+        t[i] = (v[0] * static_cast<float>(p[0]) + v[1] * static_cast<float>(p[1]) + v[2] * static_cast<float>(p[2])) * 0.5f;
+        const UnitCell got = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell();
+        const double mismatch = std::max({std::fabs(got.a / want.a - 1.0), std::fabs(got.b / want.b - 1.0),
+                                          std::fabs(got.c / want.c - 1.0),
+                                          std::fabs(got.alpha - want.alpha) / 90.0, std::fabs(got.beta - want.beta) / 90.0,
+                                          std::fabs(got.gamma - want.gamma) / 90.0});
+        if (mismatch < best_mismatch) {
+            best_mismatch = mismatch;
+            out.parity = p;
+        }
+    }
+    // A sub-lattice that is the half cell to within a percent; nothing else is the class.
+    if (best_mismatch > 0.01)
+        return out;
+
+    // Each added reflection in units of the mean of the lattice's own reflections in its shell
+    // (SUPERCELL_TREND_D_A, the shells the probe reads), in the shells where that mean is measured -
+    // three standard errors from zero - and that hold at least MIN_SHELL added reflections.
+    constexpr int64_t MIN_SHELL = 200;
+    const auto odd = [&](const MergedReflection &r) {
+        return (std::abs(out.parity[0] * r.h + out.parity[1] * r.k + out.parity[2] * r.l) & 1) != 0;
+    };
+    const auto shell_of = [](float d) {
+        if (!(d <= SUPERCELL_TREND_D_A[0]))
+            return -1;
+        int sh = 0;
+        while (d <= SUPERCELL_TREND_D_A[sh + 1])
+            sh++;
+        return sh;
+    };
+    std::array<double, SUPERCELL_TREND_SHELLS> sum{}, sum2{};
+    std::array<int64_t, SUPERCELL_TREND_SHELLS> cnt{}, cnt_odd{};
+    for (const auto &r : merged) {
+        const int sh = shell_of(r.d);
+        if (sh < 0 || !std::isfinite(r.I))
+            continue;
+        if (odd(r))
+            cnt_odd[sh]++;
+        else {
+            sum[sh] += r.I; sum2[sh] += static_cast<double>(r.I) * r.I; cnt[sh]++;
+        }
+    }
+    std::array<bool, SUPERCELL_TREND_SHELLS> use{};
+    for (int sh = 0; sh < SUPERCELL_TREND_SHELLS; sh++) {
+        if (cnt[sh] < 2 || cnt_odd[sh] < MIN_SHELL)
+            continue;
+        const double n = static_cast<double>(cnt[sh]), m = sum[sh] / n;
+        use[sh] = m > 3.0 * std::sqrt(std::max(0.0, sum2[sh] / n - m * m) / n);
+    }
+    struct Row { double h, k, l, y; int shell; };
+    std::vector<Row> rows;
+    for (const auto &r : merged) {
+        const int sh = shell_of(r.d);
+        if (sh >= 0 && use[sh] && std::isfinite(r.I) && odd(r))
+            rows.push_back({static_cast<double>(r.h), static_cast<double>(r.k), static_cast<double>(r.l),
+                            r.I / (sum[sh] / static_cast<double>(cnt[sh])), sh});
+    }
+    if (rows.size() < 1000)
+        return out;
+
+    // Least squares of y on sin^2(pi h.eps) and a constant B, for one eps, over the reflections of one
+    // shell (or all of them, shell < 0); `stride` reads every stride-th reflection only - the coarse
+    // grid needs the correlation's shape, not its last digit. Returns the correlation, which the eps
+    // search maximises, and gives the shell's mean y and B.
+    const auto fit = [&](const std::array<double, 3> &e, size_t stride, int shell, double *mean, double *b,
+                         double *n_out) {
+        double sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, n = 0;
+        for (size_t j = 0; j < rows.size(); j += stride) {
+            const auto &r = rows[j];
+            if (shell >= 0 && r.shell != shell)
+                continue;
+            n += 1;
+            const double x = std::pow(std::sin(PI * (r.h * e[0] + r.k * e[1] + r.l * e[2])), 2);
+            sx += x; sy += r.y; sxx += x * x; sxy += x * r.y; syy += r.y * r.y;
+        }
+        const double vxx = sxx - sx * sx / n, vxy = sxy - sx * sy / n, vyy = syy - sy * sy / n;
+        if (n < 3 || vxx <= 0.0 || vyy <= 0.0)
+            return -1.0;
+        if (mean) {
+            *mean = sy / n;
+            *b = (sy - vxy / vxx * sx) / n;
+            *n_out = n;
+        }
+        return vxy / std::sqrt(vxx * vyy);
+    };
+    // eps from a coarse grid over the translations within an eighth of a cell of the half point, then
+    // two finer grids about the best.
+    std::array<double, 3> best{};
+    double best_cc = -INFINITY;
+    const size_t coarse_stride = std::max<size_t>(1, rows.size() / 4000);
+    for (const auto &[step, reach] : {std::pair{0.01, 0.12}, std::pair{0.002, 0.012}, std::pair{0.0005, 0.003}}) {
+        const size_t stride = step == 0.01 ? coarse_stride : 1;
+        const std::array<double, 3> centre = best;
+        const int nstep = static_cast<int>(std::lround(reach / step));
+        for (int x = -nstep; x <= nstep; x++)
+            for (int y = -nstep; y <= nstep; y++)
+                for (int z = -nstep; z <= nstep; z++) {
+                    const std::array<double, 3> e = {centre[0] + x * step, centre[1] + y * step, centre[2] + z * step};
+                    const double cc = fit(e, stride, -1, nullptr, nullptr, nullptr);
+                    if (cc > best_cc) {
+                        best_cc = cc;
+                        best = e;
+                    }
+                }
+    }
+    // The share of each shell's added class that the translation modulates: 1 - B over the class's mean.
+    for (int sh = 0; sh < SUPERCELL_TREND_SHELLS; sh++) {
+        if (!use[sh])
+            continue;
+        double mean = 0, b = 0, n = 0;
+        if (fit(best, 1, sh, &mean, &b, &n) < -0.5 || !(mean > 0.0))
+            continue;
+        out.shells.push_back({SUPERCELL_TREND_D_A[sh], SUPERCELL_TREND_D_A[sh + 1], static_cast<int64_t>(n),
+                              mean, 1.0 - b / mean});
+    }
+    out.eps = best;
+    out.tested = !out.shells.empty();
+    return out;
+}
+
+namespace {
+    std::string CellText(const UnitCell &c) {
+        return fmt::format("{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}", c.a, c.b, c.c, c.alpha, c.beta, c.gamma);
+    }
+}
+
+ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessResult half,
+                                            const std::string &base_prefix) {
+    if (cancelled_ || user_cell_ || user_fixed_sg_ || !half.supercell || !last_rotation_result_
+        || !half.has_merge_statistics || !half.consensus_cell || !DoubledCellWorthMerging(*half.supercell))
+        return half;
+    Logger logger("Rugnux");
+    const auto &sc = *half.supercell;
+    // Only a class that could be adopted is worth the pass: one that does not rise with resolution keeps
+    // the half cell whatever the doubled merge shows, and a doubled merge of a large crystal is not free
+    // (measured: an R-centred 531 A cell ran the card out of memory).
+    const std::string rises = DoubledCellRises(sc);
+    if (rises.empty()) {
+        logger.Info("Doubled-cell hypothesis: the half-integer class {} {} {} stands clear of the other classes but "
+                    "does not grow with resolution to a quarter of the lattice's intensity, so the half cell describes "
+                    "the crystal at high resolution too - not merged", sc.h, sc.k, sc.l);
+        return half;
+    }
+
+    // The lattice the probed cell was integrated on, doubled along the class, with everything else of
+    // that indexing - orientation, geometry, axis - kept, so the two merges differ by the cell alone.
+    RotationIndexerResult doubled = *last_rotation_result_;
+    const CrystalLattice prim = doubled.lattice.ToPrimitive(doubled.search_result.centering);
+    // Classed only on what the refined metric already holds exactly: a looser search promotes a doubled
+    // cell whose angle merely lies near 90 deg (measured: 88.2 deg, held to 90, and no frame indexed).
+    doubled.search_result = LatticeSearch(DoubledLattice(prim, {sc.h, sc.k, sc.l}), 0.001, 0.05);
+    doubled.lattice = doubled.search_result.conventional;
+    doubled.extra_lattices.clear();
+    doubled.unconstrained.reset();
+    logger.Info("Doubled-cell hypothesis: the half-integer class {} {} {} stands clear of the other classes out "
+                "to the limit, so the run integrates and merges the doubled cell as well - lattice {} ({}-centred "
+                "{}) doubled to {} ({}-centred {})", sc.h, sc.k, sc.l,
+                CellText(last_rotation_result_->lattice.GetUnitCell()), last_rotation_result_->search_result.centering,
+                gemmi::crystal_system_str(last_rotation_result_->search_result.system),
+                CellText(doubled.lattice.GetUnitCell()), doubled.search_result.centering,
+                gemmi::crystal_system_str(doubled.search_result.system));
+    if (observer)
+        observer->OnStep("Doubled-cell hypothesis");
+
+    const std::string half_prefix = base_prefix + "_halfcell";
+    const std::string doubled_prefix = base_prefix + "_doubled";
+    config_.output_prefix = base_prefix.empty() ? base_prefix : doubled_prefix;
+    force_rotation_result_ = doubled;
+    keep_merged_reflections_ = true;
+
+    // The quality guard compares a refined pass against the header pass in the same cell; this pass is
+    // in another cell and is judged here instead.
+    quality_guard_pass1_.reset();
+    ProcessResult d;
+    try {
+        d = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
+    } catch (const std::exception &e) {
+        if (IsFatalResourceError(e)) throw;
+        logger.Warning("Doubled-cell hypothesis: the pass in the doubled cell did not complete ({})", e.what());
+    }
+    force_rotation_result_.reset();
+    keep_merged_reflections_ = false;
+    config_.output_prefix = base_prefix;
+    if (!d.has_merge_statistics || !d.consensus_cell) {
+        logger.Warning("Doubled-cell hypothesis: the doubled cell merged nothing - the probed cell stays the answer");
+        return half;
+    }
+    d.pass_number = d.pass_count = half.pass_count + 1;
+    half.pass_count = d.pass_count;
+
+    // Decided with both merges in hand. The doubled cell is the answer where the added class is a
+    // pseudo-translation breaking down - it grows with resolution (DoubledCellRises) - and where, in the
+    // doubled merge, that one translation accounts for most of the class in every shell: the two halves
+    // of the cell are then one content displaced, out to the limit (FitDoubledClassTranslation).
+    const DoubledClassTranslation tr = FitDoubledClassTranslation(d.merged_reflections, *d.consensus_cell,
+                                                                  *half.consensus_cell,
+                                                                  half.consensus_centering.value_or('P'));
+    d.merged_reflections.clear();
+    const DoubledClassTranslation::Shell *weakest = nullptr;
+    std::string table;
+    for (const auto &sh : tr.shells) {
+        table += fmt::format(" {:.1f}-{:.1f} A {:.2f};", sh.d_low_A, sh.d_high_A, sh.modulated_share);
+        if (!weakest || sh.modulated_share < weakest->modulated_share)
+            weakest = &sh;
+    }
+    if (tr.tested)
+        logger.Info("Doubled-cell hypothesis: the added class {} {} {} of the doubled merge against one pseudo-"
+                    "translation, half a cell plus ({:.4f}, {:.4f}, {:.4f}) - the share of it the translation "
+                    "modulates, by shell:{}", tr.parity[0], tr.parity[1], tr.parity[2], tr.eps[0], tr.eps[1],
+                    tr.eps[2], table);
+    const bool one_translation = tr.tested && weakest->modulated_share >= 0.5;
+    // A superstructure cannot be more symmetric than its own average: the half cell's merge is that
+    // average, so a doubled merge that reads a point group of higher order than the half cell's describes
+    // something other than the crystal's superstructure (measured: a C2221 half cell whose doubled merge
+    // came out centrosymmetric, Cccm, with the added class otherwise passing every test).
+    const auto point_group_order = [](const ProcessResult &r) {
+        return r.space_group ? r.space_group->operations().sym_ops.size() : size_t{0};
+    };
+    const bool symmetry_raised = point_group_order(d) > point_group_order(half);
+    const bool adopt = one_translation && !symmetry_raised;
+    ProcessResult::DoubledCell h;
+    h.adopted = adopt;
+    h.decision = adopt
+            ? fmt::format("{}; and one pseudo-translation accounts for most of it in every shell (at least {:.2f}, "
+                          "at {:.1f}-{:.1f} A)", rises, weakest->modulated_share, weakest->d_low_A, weakest->d_high_A)
+            : symmetry_raised
+            ? fmt::format("the probed cell stands: the doubled merge reads {}, more symmetric than the half cell's {}, "
+                          "which a superstructure of it cannot be", d.space_group->xhm(), half.space_group->xhm())
+            : !tr.tested
+            ? std::string("the probed cell stands: the added class could not be read in the doubled merge")
+            : fmt::format("the probed cell stands: one pseudo-translation accounts for only {:.2f} of the added class "
+                          "at {:.1f}-{:.1f} A - the rest is there whatever the direction, which copies related by a "
+                          "translation do not put there", weakest->modulated_share, weakest->d_low_A,
+                          weakest->d_high_A);
+    logger.Info("Doubled-cell hypothesis: {}", h.decision);
+
+    // The answer's files go under the plain prefix and the other cell's beside it: where the doubled cell
+    // is adopted, every file it wrote takes the plain name, and the probed cell's file of the same name
+    // moves to <prefix>_halfcell.
+    if (adopt && !base_prefix.empty()) {
+        const std::filesystem::path dp(doubled_prefix);
+        const std::filesystem::path dir = dp.has_parent_path() ? dp.parent_path() : std::filesystem::path(".");
+        const std::string stem = dp.filename().string();
+        std::error_code ec;
+        std::vector<std::string> suffixes;
+        for (const auto &f : std::filesystem::directory_iterator(dir, ec))
+            if (const std::string name = f.path().filename().string(); name.starts_with(stem))
+                suffixes.push_back(name.substr(stem.size()));
+        for (const auto &x : suffixes) {
+            if (std::filesystem::exists(base_prefix + x, ec))
+                std::filesystem::rename(base_prefix + x, half_prefix + x, ec);
+            std::filesystem::rename(doubled_prefix + x, base_prefix + x, ec);
+        }
+    }
+
+    ProcessResult &answer = adopt ? d : half;
+    const ProcessResult &other = adopt ? half : d;
+    h.other_prefix = adopt ? half_prefix : doubled_prefix;
+    h.other_cell = *other.consensus_cell;
+    h.other_space_group = other.space_group ? other.space_group->xhm() : std::string();
+    h.other_d_min_A = other.merge_statistics.overall.d_min;
+    h.other_cc_half = other.merge_statistics.overall.cc_half;
+    if (other.model_validation)
+        h.other_r_free = other.model_validation->r_free;
+    answer.doubled_cell = h;
+    // The answer's report describes the probe of the probed cell, which is what the decision read.
+    answer.supercell = half.supercell;
+    answer.pass_decision = fmt::format("doubled-cell hypothesis merged ({}); {}", h.decision, half.pass_decision);
+    return answer;
+}

rugnux/RugnuxScaleMerge.cpp

@@ -3983,6 +3983,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
         WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics,
                          em_report, result.twinning, config_.provenance, config_.output_prefix,
                          static_cast<size_t>(std::max(1, config_.nthreads)), &sm.scaled_fulls);
+        if (keep_merged_reflections_)
+            result.merged_reflections = sm.merged;
         // The twin hypothesis beside it (see prepare_twin_file): the reflection files only, with the
         // law and its fraction in the mmCIF's _pdbx_reflns_twin.
         if (twin_file) {

tests/IndexingUnitTest.cpp

@@ -817,6 +817,23 @@ TEST_CASE("FitSupercellProbe", "[SupercellProbe]") {
     CHECK(FitSupercellProbe(SupercellProbeClass{}).b == 0.0);
 }
 
+TEST_CASE("DoubledLattice", "[SupercellProbe]") {
+    const CrystalLattice prim(Coord(40, 0, 0), Coord(0, 50, 0), Coord(0, 0, 60));
+    // Class 0 0 1: the c axis doubled, the others kept.
+    const CrystalLattice c2 = DoubledLattice(prim, {0, 0, 1});
+    CHECK(c2.Vec0().x == Catch::Approx(40.0));
+    CHECK(c2.Vec1().y == Catch::Approx(50.0));
+    CHECK(c2.Vec2().z == Catch::Approx(120.0));
+    // Class 1 1 0: a centred doubling of the ab plane - twice a, and b plus a - so every vector of the
+    // supercell has h + k even, at twice the volume.
+    const CrystalLattice ab = DoubledLattice(prim, {1, 1, 0});
+    CHECK(std::fabs(ab.CalcVolume()) == Catch::Approx(2.0 * 40 * 50 * 60));
+    CHECK(ab.Vec0().x == Catch::Approx(80.0));
+    CHECK(ab.Vec1().x == Catch::Approx(40.0));
+    CHECK(ab.Vec1().y == Catch::Approx(50.0));
+    CHECK(ab.Vec2().z == Catch::Approx(60.0));
+}
+
 TEST_CASE("PostIndexingRefinement_CoplanarCandidateIsRejected", "[Indexing]") {
     // A candidate whose three rows lie in one plane passes every length and angle test - three rows
     // 120 deg apart in a plane meet the 30-150 deg window - and every spot along the plane normal

tests/RugnuxTest.cpp

@@ -1,6 +1,7 @@
 // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only
 
+#include <random>
 #include <catch2/catch_all.hpp>
 #include <cmath>
 #include <filesystem>
@@ -349,6 +350,83 @@ TEST_CASE("SpotWidth_Gaussian", "[process]") {
 // The rotation two-pass quality guard. The refined pass is sent back to the header geometry only for
 // a decisive loss of signal - reflections at I/sigma >= 2 net of the noise tail at <= -2 - over no
 // more of reciprocal space.
+TEST_CASE("DoubledCellDecision", "[process]") {
+    // Shells of the supercell probe's trend: d range, excess over the null, z, the strongest null z and
+    // the lattice's <I/sigma>.
+    const auto trend = [](std::vector<std::array<double, 3>> excess_z_lattice) {
+        ProcessResult::SupercellClass sc;
+        float d = 20.0f;
+        for (const auto &[excess, z, lattice] : excess_z_lattice) {
+            ProcessResult::SupercellClass::Shell sh;
+            sh.d_low_A = d;
+            sh.d_high_A = d *= 0.7f;
+            sh.excess = excess;
+            sh.z = z;
+            sh.sibling_z = 3.0;
+            sh.lattice_i_over_sigma = lattice;
+            sh.resolved = 0.9;
+            sc.trend.push_back(sh);
+        }
+        return sc;
+    };
+    // A pseudo-translation that breaks down: weak at low resolution, a third of the lattice where the
+    // lattice is still measured.
+    const auto breaks_down = trend({{0.02, 17, 18}, {0.10, 38, 7}, {0.26, 36, 6}, {0.35, 45, 2.4}, {0.46, 33, 0.8}});
+    CHECK(DoubledCellWorthMerging(breaks_down));
+    CHECK(!DoubledCellRises(breaks_down).empty());
+    // A class that stays a tenth of the lattice at every resolution: merged as well, but it does not rise.
+    const auto flat = trend({{0.09, 10, 70}, {0.11, 23, 28}, {0.10, 24, 33}, {0.11, 24, 18}, {0.10, 23, 8}});
+    CHECK(DoubledCellWorthMerging(flat));
+    CHECK(DoubledCellRises(flat).empty());
+    // A class no stronger than the null: nothing.
+    const auto null_class = trend({{0.01, 4, 20}, {0.02, 5, 10}, {0.03, 2, 5}, {0.05, 1, 2}});
+    CHECK(!DoubledCellWorthMerging(null_class));
+    // A class that rises significantly but stays a few percent of the lattice: the half cell holds.
+    const auto small = trend({{0.009, 40, 500}, {0.014, 60, 300}, {0.016, 90, 170}, {0.03, 60, 50}, {0.06, 40, 2}});
+    CHECK(DoubledCellRises(small).empty());
+}
+
+TEST_CASE("FitDoubledClassTranslation", "[process]") {
+    // A doubled orthorhombic cell whose two halves are one copy translated by c/2 + (0.03, 0, 0): the
+    // reflections with l odd carry |F|^2 sin^2(pi h.eps), the others cos^2. With nothing else on the
+    // added class the fit finds the translation and leaves nothing over; with a flat 0.3 added to
+    // the class, most of the class is left unexplained.
+    const UnitCell doubled{50, 60, 140, 90, 90, 90}, half{50, 60, 70, 90, 90, 90};
+    std::mt19937 rng(1);
+    std::exponential_distribution<double> wilson(1.0);
+    std::vector<MergedReflection> clean, contaminated;
+    for (int h = -25; h <= 25; h++)
+        for (int k = 0; k <= 30; k++)
+            for (int l = 0; l <= 70; l++) {
+                if (h == 0 && k == 0 && l == 0)
+                    continue;
+                MergedReflection r;
+                r.h = h; r.k = k; r.l = l;
+                r.d = 1.0f / std::sqrt(static_cast<float>(h * h / 2500.0 + k * k / 3600.0 + l * l / 19600.0));
+                if (r.d < 2.0f)
+                    continue;
+                const double base = wilson(rng), x = std::pow(std::sin(PI * 0.03 * h), 2);
+                r.I = static_cast<float>(base * ((l & 1) ? x : 1.0 - x));
+                clean.push_back(r);
+                if (l & 1)
+                    r.I += 0.3f;
+                contaminated.push_back(r);
+            }
+    const DoubledClassTranslation a = FitDoubledClassTranslation(clean, doubled, half);
+    REQUIRE(a.tested);
+    CHECK(a.parity == std::array<int, 3>{0, 0, 1});
+    const DoubledClassTranslation b = FitDoubledClassTranslation(contaminated, doubled, half);
+    REQUIRE(b.tested);
+    const auto weakest = [](const DoubledClassTranslation &t) {
+        double w = INFINITY;
+        for (const auto &sh : t.shells)
+            w = std::min(w, sh.modulated_share);
+        return w;
+    };
+    CHECK(weakest(a) > 0.8);
+    CHECK(weakest(b) < 0.5);
+}
+
 TEST_CASE("RefinedPassIsWorse", "[process]") {
     const auto pass = [](const UnitCell &cell, char centering, int64_t strong, int64_t negative,
                          int64_t reflections) {