WIP: rugnux doubled cell decided late on one pseudo-translation; screw trim reads the measured intensity

Doubled-cell hypothesis: the run now always merges the doubled cell under <prefix>_doubled first and
decides with both merges in hand. Adopted only where (1) the probe's class grows significantly with
resolution (three standard errors between the 20-8 A shell and the finest measured one) and (2) in
the doubled merge one pseudo-translation - half a cell plus eps, fitted over the class - modulates at
least half of the added class in every shell (FitDoubledClassTranslation: A sin^2(pi h.eps) + B per
shell, share 1 - B/mean). Where adopted, every <prefix>_doubled file takes the plain name and the
probed cell's file of that name moves to <prefix>_halfcell. Replaces the HALF_CELL_BREAKS_DOWN bound.

Measured on the forced doubled merges: modulated share per shell 0.80-1.0 out to 1.9 A on the crystal
deposited in the doubled cell; on the ones deposited in the half cell it falls to 0.01-0.34 by
2.1-2.8 A (the rising class with twinning), or stays at -0.3-0.48 at every resolution (three others);
the one class it fully explains does not rise (eps 0.1 of a cell, flat at 0.1 of the lattice).

Space-group search: the trim of a screw zone's largest absence exempted a violation "at its row's
mean" on the de-modulated intensity, so under a deep pseudo-translation a reflection at 0.015 of its
row (de-modulated by 1/0.002) could never be trimmed and alone decided the row against the screw.
The exemption now reads the measured intensity, at both sites.

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

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/Rugnux.h

@@ -530,6 +530,9 @@ struct ProcessResult {
         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
@@ -650,10 +653,32 @@ std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &
 
 // 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.
-// Adopt: the doubled cell is the answer, under DOUBLED_CELL_POLICY; empty when the probed cell stands,
-// otherwise why the doubled one is adopted.
+// 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 AdoptDoubledCell(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);
 
 // Whether two passes hold the same lattice: the same Bravais class and primitive volumes within 2 %.
 bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b);
@@ -769,6 +794,8 @@ class Rugnux {
     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,

rugnux/RugnuxPasses.cpp

@@ -1415,27 +1415,24 @@ std::optional<ProcessResult> Rugnux::RefinedPasses(RugnuxObserver *observer, Pro
 //
 // 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),
-// but nothing read from intensities told the crystal whose deposited cell is the doubled one from those
-// deposited in the half cell - the odd/even ratio rose with resolution alike on both kinds, and the
-// odd class correlated with its lattice neighbours alike. Which of the two cells to write is therefore
-// a policy (DOUBLED_CELL_POLICY), and both are written.
+// 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.
     bool DoubledCellShellMeasured(const ProcessResult::SupercellClass::Shell &sh) {
         return sh.z >= 3.0 && sh.lattice_i_over_sigma >= 1.0;
     }
-
-    // ADDED_CLASS_SIGNIFICANT: the doubled cell wherever the class stands clear of the null - it
-    //   describes every Bragg reflection the crystal has, where the half cell leaves a measured class
-    //   out.
-    // HALF_CELL_BREAKS_DOWN: the doubled cell only where the class reaches DOUBLED_CELL_BREAKDOWN_EXCESS
-    //   of the lattice's intensity in the finest shell both are measured in - where the half cell no
-    //   longer approximates the crystal. A class that stays a small fraction of the lattice at every
-    //   resolution (two halves that differ slightly throughout) keeps the half cell, as the depositions
-    //   of such crystals do.
-    enum class DoubledCellPolicy { AddedClassSignificant, HalfCellBreaksDown };
-    constexpr DoubledCellPolicy DOUBLED_CELL_POLICY = DoubledCellPolicy::HalfCellBreaksDown;
-    constexpr double DOUBLED_CELL_BREAKDOWN_EXCESS = 0.25;
 }
 
 bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc) {
@@ -1447,18 +1444,157 @@ bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc) {
     return clear >= 3;
 }
 
-std::string AdoptDoubledCell(const ProcessResult::SupercellClass &sc) {
-    if (DOUBLED_CELL_POLICY == DoubledCellPolicy::AddedClassSignificant)
-        return DoubledCellWorthMerging(sc) ? "the added class stands clear of the other classes" : std::string();
+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: a significant rise, not a size.
+    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;
-    if (!outer || outer->excess < DOUBLED_CELL_BREAKDOWN_EXCESS)
+    const double se = std::hypot(low.excess / low.z, outer->excess / outer->z);
+    if (outer == &low || outer->excess - low.excess < 3.0 * se)
         return {};
-    return fmt::format("the added class reaches {:.2f} of the lattice's intensity at {:.1f}-{:.1f} A, the finest "
-                       "shell both are measured in, so the half cell no longer approximates the crystal",
-                       outer->excess, outer->d_low_A, outer->d_high_A);
+    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) {
+    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).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 {
@@ -1474,8 +1610,6 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes
         return half;
     Logger logger("Rugnux");
     const auto &sc = *half.supercell;
-    const std::string adopt_why = AdoptDoubledCell(sc);
-    const bool adopt = !adopt_why.empty();
 
     // 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.
@@ -1487,32 +1621,21 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes
     doubled.lattice = doubled.search_result.conventional;
     doubled.extra_lattices.clear();
     doubled.unconstrained.reset();
-    logger.Info("Doubled-cell hypothesis: lattice {} ({}-centred {}) doubled to {} ({}-centred {})",
+    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));
-
-    // The answer's files go under the plain prefix and the other cell's merged files beside them, so
-    // where the doubled cell is adopted the probed cell's merged files move aside before it writes.
-    const std::string half_prefix = base_prefix + "_halfcell";
-    const std::string doubled_prefix = base_prefix + "_doubled";
-    const auto move_merged = [&](const std::string &from, const std::string &to) {
-        for (const char *ext : {".mtz", ".cif", ".hkl"})
-            if (std::error_code ec; std::filesystem::exists(from + ext, ec))
-                std::filesystem::rename(from + ext, to + ext, ec);
-    };
-    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 {:.2f} {:.2f} {:.2f} "
-                "{:.2f} {:.2f} {:.2f} as well; {}", sc.h, sc.k, sc.l, sc.doubled_cell.a, sc.doubled_cell.b,
-                sc.doubled_cell.c, sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma,
-                adopt ? "it is the answer - " + adopt_why : std::string("the probed cell stays the answer"));
     if (observer)
         observer->OnStep("Doubled-cell hypothesis");
-    if (adopt && !base_prefix.empty())
-        move_merged(base_prefix, half_prefix);
-    config_.output_prefix = base_prefix.empty() ? base_prefix : adopt ? base_prefix : doubled_prefix;
+
+    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();
@@ -1524,24 +1647,74 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes
         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");
-        if (adopt && !base_prefix.empty())
-            move_merged(half_prefix, base_prefix);
         return half;
     }
     d.pass_number = d.pass_count = half.pass_count + 1;
     half.pass_count = d.pass_count;
 
-    ProcessResult &answer = adopt ? d : half;
-    const ProcessResult &other = adopt ? half : d;
+    // 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 std::string rises = DoubledCellRises(sc);
+    const DoubledClassTranslation tr = FitDoubledClassTranslation(d.merged_reflections, *d.consensus_cell,
+                                                                  *half.consensus_cell);
+    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;
+    const bool adopt = !rises.empty() && one_translation;
     ProcessResult::DoubledCell h;
     h.adopted = adopt;
-    h.decision = adopt ? adopt_why : DOUBLED_CELL_POLICY == DoubledCellPolicy::HalfCellBreaksDown
-            ? fmt::format("the probed cell stands: the added class stays below {:.2f} of the lattice's intensity "
-                          "where both are measured", DOUBLED_CELL_BREAKDOWN_EXCESS)
-            : std::string("the probed cell stands");
+    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)
+            : rises.empty()
+            ? std::string("the probed cell stands: the added class does not grow significantly with resolution, "
+                          "so the half cell describes the crystal as well at high resolution as at low")
+            : !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();

rugnux/RugnuxScaleMerge.cpp

@@ -3974,6 +3974,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/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>
@@ -371,16 +372,57 @@ TEST_CASE("DoubledCellDecision", "[process]") {
     // 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(!AdoptDoubledCell(breaks_down).empty());
-    // A class that stays a tenth of the lattice at every resolution: merged as well, the half cell kept.
+    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(AdoptDoubledCell(flat).empty());
+    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));
 }
 
+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) {