Doubled-cell hypothesis also on sweeps whose probe frames do not measure the lattice

The doubled-cell hypothesis needed the supercell probe to measure the lattice's own reflections
(per-observation <I/sigma> of 1 or more) in the shells where the half-integer class stands clear of
the six empty classes. On a weak sweep no probe shell measures the lattice, so a class that clears
its siblings was never followed up. Such a run now merges the doubled cell as well where the class
clears the siblings in at least two shells, and reads the class's rise on the doubled merge (the
added class against the lattice's own mean intensity per probe shell) with the same rise test; the
translation and symmetry tests then decide as before. Runs where any probe shell measures the
lattice are unchanged. The report says where the rise was read (DOUBLED_CELL_RISE_READ_ON).

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

docs/RUGNUX_REPORT.md

@@ -739,6 +739,11 @@ small offset, fitted over the class - modulates at least half of the added class
 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.
+On a sweep too weak for the probe's 60 frames to measure the lattice's own reflections in any shell, the
+probe cannot read the rise; there the class only has to stand clear of the six empty ones in two shells,
+the doubled cell is merged, and the rise is read on the doubled merge instead - the added class against
+the lattice's own mean intensity, shell by shell, with the same test - before the same translation and
+symmetry tests decide. **`DOUBLED_CELL_RISE_READ_ON=`** says which (`PROBE` or `MERGE`).
 **`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

rugnux/ResultReport.cpp

@@ -1847,6 +1847,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
             if (result.doubled_cell) {
                 const auto &h = *result.doubled_cell;
                 Add(s, KeyBool("DOUBLED_CELL_ADOPTED", h.adopted));
+                Add(s, KeyEnum("DOUBLED_CELL_RISE_READ_ON", h.rise_on_merge ? "MERGE" : "PROBE", {"PROBE", "MERGE"}));
                 Add(s, KeyText("DOUBLED_CELL_DECISION", h.decision));
                 Add(s, KeyText("DOUBLED_CELL_OTHER_FILE", output_prefix + h.other_suffix + ".mtz"));
                 Add(s, KeyText("DOUBLED_CELL_OTHER_CELL", CellString(h.other_cell)));

rugnux/Rugnux.h

@@ -768,6 +768,7 @@ struct ProcessResult {
     // 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
+        bool rise_on_merge = false;        // the class's rise was read on the doubled merge, not the probe
         std::string decision;              // why, in the words the run reports it in
         std::string other_suffix;          // _halfcell or _doubled: the other cell's files are <prefix><other_suffix>
         std::shared_ptr<RugnuxOutputs> other_outputs;
@@ -974,6 +975,14 @@ std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &
 // 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);
+// The same hypothesis on a sweep too weak for the probe's frames to measure the lattice in any shell: worth
+// merging where the class stands clear of the six empty ones in at least two shells, and the rise is then
+// read on the doubled merge - the added class against the lattice's own, per probe shell.
+bool DoubledCellWorthMergingWeak(const ProcessResult::SupercellClass &sc);
+std::optional<std::array<int, 3>> DoubledClassParity(const UnitCell &doubled_cell, const UnitCell &half_cell,
+                                                     char half_centering);
+std::vector<ProcessResult::SupercellClass::Shell> MergedDoubledClassTrend(const std::vector<MergedReflection> &merged,
+                                                                          const std::array<int, 3> &parity);
 
 // 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

rugnux/RugnuxPasses.cpp

@@ -1695,10 +1695,23 @@ std::string DoubledCellRises(const ProcessResult::SupercellClass &sc) {
                        (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;
+bool DoubledCellWorthMergingWeak(const ProcessResult::SupercellClass &sc) {
+    // On a weak sweep the probe's frames do not measure the lattice's own reflections in any shell (a single
+    // observation reads <I/sigma> below 1), so DoubledCellWorthMerging cannot speak whatever the class does.
+    // The class is still read against the six empty ones, pooled over the same frames: where it stands
+    // clear of them in at least two shells, the doubled merge is made and the class is read there instead
+    // (MergedDoubledClassTrend).
+    int clear = 0;
+    for (const auto &sh : sc.trend) {
+        if (sh.lattice_i_over_sigma >= 1.0)
+            return false;
+        clear += sh.z >= 3.0 && sh.z > 2.0 * sh.sibling_z && sh.resolved >= 0.5;
+    }
+    return clear >= 2;
+}
+
+std::optional<std::array<int, 3>> DoubledClassParity(const UnitCell &doubled_cell, const UnitCell &half_cell,
+                                                     char half_centering) {
     // 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.
@@ -1706,6 +1719,7 @@ DoubledClassTranslation FitDoubledClassTranslation(const std::vector<MergedRefle
     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;
+    std::array<int, 3> parity{};
     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;
@@ -1718,12 +1732,63 @@ DoubledClassTranslation FitDoubledClassTranslation(const std::vector<MergedRefle
                                           std::fabs(got.gamma - want.gamma) / 90.0});
         if (mismatch < best_mismatch) {
             best_mismatch = mismatch;
-            out.parity = p;
+            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 std::nullopt;
+    return parity;
+}
+
+std::vector<ProcessResult::SupercellClass::Shell> MergedDoubledClassTrend(const std::vector<MergedReflection> &merged,
+                                                                          const std::array<int, 3> &parity) {
+    // The probe's trend (SupercellClass::trend) read on the doubled merge: per shell, the added class's
+    // mean intensity against the lattice's own (excess) and in its own standard errors (z), and the
+    // lattice's merged <I/sigma>. Every added reflection of the merge is integrated in its own right, so
+    // all of the class is resolved; a doubled cell has no sibling classes to set beside it.
+    std::vector<ProcessResult::SupercellClass::Shell> out;
+    std::array<double, SUPERCELL_TREND_SHELLS> s_odd{}, s2_odd{}, s_even{}, s_ios{};
+    std::array<int64_t, SUPERCELL_TREND_SHELLS> n_odd{}, n_even{};
+    for (const auto &r : merged) {
+        if (!(r.d <= SUPERCELL_TREND_D_A[0]) || !std::isfinite(r.I) || !(r.sigma > 0.0f))
+            continue;
+        int sh = 0;
+        while (sh + 1 < SUPERCELL_TREND_SHELLS && r.d <= SUPERCELL_TREND_D_A[sh + 1])
+            sh++;
+        if ((std::abs(parity[0] * r.h + parity[1] * r.k + parity[2] * r.l) & 1) != 0) {
+            s_odd[sh] += r.I; s2_odd[sh] += static_cast<double>(r.I) * r.I; n_odd[sh]++;
+        } else {
+            s_even[sh] += r.I; s_ios[sh] += r.I / r.sigma; n_even[sh]++;
+        }
+    }
+    for (int sh = 0; sh < SUPERCELL_TREND_SHELLS; sh++) {
+        if (n_odd[sh] < 2 || n_even[sh] < 2 || !(s_even[sh] > 0.0))
+            continue;
+        const double no = static_cast<double>(n_odd[sh]), ne = static_cast<double>(n_even[sh]);
+        const double mean_even = s_even[sh] / ne, mean_odd = s_odd[sh] / no;
+        const double se = std::sqrt(std::max(0.0, s2_odd[sh] / no - mean_odd * mean_odd) / no) / mean_even;
+        ProcessResult::SupercellClass::Shell row;
+        row.d_low_A = SUPERCELL_TREND_D_A[sh];
+        row.d_high_A = SUPERCELL_TREND_D_A[sh + 1];
+        row.n = n_odd[sh];
+        row.excess = mean_odd / mean_even;
+        row.z = se > 0.0 ? row.excess / se : 0.0;
+        row.lattice_i_over_sigma = s_ios[sh] / ne;
+        row.resolved = 1.0;
+        out.push_back(row);
+    }
+    return out;
+}
+
+DoubledClassTranslation FitDoubledClassTranslation(const std::vector<MergedReflection> &merged,
+                                                   const UnitCell &doubled_cell, const UnitCell &half_cell,
+                                                   char half_centering) {
+    DoubledClassTranslation out;
+    const auto parity = DoubledClassParity(doubled_cell, half_cell, half_centering);
+    if (!parity)
         return out;
+    out.parity = *parity;
 
     // 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 -
@@ -1838,15 +1903,23 @@ namespace {
 
 ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessResult half) {
     if (cancelled_ || config_.fast_live || user_cell_ || user_fixed_sg_ || !half.supercell || !last_rotation_result_
-        || !half.has_merge_statistics || !half.consensus_cell || !DoubledCellWorthMerging(*half.supercell))
+        || !half.has_merge_statistics || !half.consensus_cell)
+        return half;
+    // Two ways in. The probe measures the class against a lattice it measures too, and the rise is read
+    // there (DoubledCellWorthMerging); or the sweep is too weak for the probe's frames to measure the
+    // lattice at all, and the class only stands clear of the six empty ones - then the rise is read on the
+    // doubled merge (DoubledCellWorthMergingWeak, MergedDoubledClassTrend).
+    const bool on_probe = DoubledCellWorthMerging(*half.supercell);
+    const bool on_merge = !on_probe && DoubledCellWorthMergingWeak(*half.supercell);
+    if (!on_probe && !on_merge)
         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()) {
+    std::string rises = on_probe ? DoubledCellRises(sc) : std::string();
+    if (on_probe && 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);
@@ -1904,6 +1977,21 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes
     const DoubledClassTranslation tr = FitDoubledClassTranslation(d.merged_reflections, *d.consensus_cell,
                                                                   *half.consensus_cell,
                                                                   half.consensus_centering.value_or('P'));
+    if (on_merge) {
+        // The rise, read on the doubled merge as the probe would have read it.
+        ProcessResult::SupercellClass merged_class = sc;
+        merged_class.trend.clear();
+        if (const auto parity = DoubledClassParity(*d.consensus_cell, *half.consensus_cell,
+                                                   half.consensus_centering.value_or('P')))
+            merged_class.trend = MergedDoubledClassTrend(d.merged_reflections, *parity);
+        std::string trend;
+        for (const auto &sh : merged_class.trend)
+            trend += fmt::format(" {:.1f}-{:.1f} A {:.3f} ({:.1f} SE, lattice <I/sigma> {:.1f});", sh.d_low_A,
+                                 sh.d_high_A, sh.excess, sh.z, sh.lattice_i_over_sigma);
+        logger.Info("Doubled-cell hypothesis: the added class read on the doubled merge, against the lattice's own "
+                    "mean intensity, by shell:{}", trend);
+        rises = DoubledCellRises(merged_class);
+    }
     d.merged_reflections.clear();
     const DoubledClassTranslation::Shell *weakest = nullptr;
     std::string table;
@@ -1926,10 +2014,14 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes
         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;
+    const bool adopt = !rises.empty() && one_translation && !symmetry_raised;
     ProcessResult::DoubledCell h;
     h.adopted = adopt;
-    h.decision = adopt
+    h.rise_on_merge = on_merge;
+    h.decision = rises.empty()
+            ? std::string("the probed cell stands: the added class, read on the doubled merge, does not grow with "
+                          "resolution to a quarter of the lattice's intensity")
+            : 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

tests/RugnuxTest.cpp

@@ -526,6 +526,49 @@ TEST_CASE("DoubledCellDecision", "[process]") {
     // 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());
+    // The weak-sweep arm speaks only where the probe measures the lattice in no shell.
+    CHECK(!DoubledCellWorthMergingWeak(breaks_down));
+    CHECK(!DoubledCellWorthMergingWeak(flat));
+    // Lattice unmeasured on the probe's frames: the class clears the null (z > 2 x 3) in two shells.
+    const auto weak = trend({{0.12, 2.4, 0.3}, {0.37, 2.8, 0.0}, {0.62, 6.8, 0.0}, {0.56, 7.0, 0.0}});
+    CHECK(!DoubledCellWorthMerging(weak));
+    CHECK(DoubledCellWorthMergingWeak(weak));
+    // ... and in one shell only: not enough.
+    const auto weak_one = trend({{0.12, 2.4, 0.3}, {0.37, 2.8, 0.0}, {0.62, 6.8, 0.0}, {0.56, 4.0, 0.0}});
+    CHECK(!DoubledCellWorthMergingWeak(weak_one));
+}
+
+TEST_CASE("MergedDoubledClassTrend", "[process]") {
+    // A doubled cell whose l-odd class is a fifth of the l-even intensity at low resolution and all of it
+    // at high: the trend reads the ratio per probe shell, and rises.
+    const UnitCell doubled{50, 60, 140, 90, 90, 90}, half{50, 60, 70, 90, 90, 90};
+    const auto parity = DoubledClassParity(doubled, half, 'P');
+    REQUIRE(parity.has_value());
+    CHECK((*parity == std::array<int, 3>{0, 0, 1}));
+    std::mt19937 rng(7);
+    std::exponential_distribution<double> wilson(1.0);
+    std::vector<MergedReflection> merged;
+    for (int h = 0; h <= 20; h++)
+        for (int k = 0; k <= 24; k++)
+            for (int l = 0; l <= 56; l++) {
+                const double s2 = h * h / 2500.0 + k * k / 3600.0 + l * l / 19600.0;
+                if (s2 <= 0.0 || s2 > 1.0 / (2.5 * 2.5))
+                    continue;
+                MergedReflection r;
+                r.h = h; r.k = k; r.l = l;
+                r.d = static_cast<float>(1.0 / std::sqrt(s2));
+                const double t = std::clamp((1.0 / 8.0 - 1.0 / r.d) / (1.0 / 8.0 - 1.0 / 2.5), 0.0, 1.0);
+                const double share = (l % 2) ? 0.2 + 0.8 * t : 1.0;
+                r.I = static_cast<float>(100.0 * share * wilson(rng));
+                r.sigma = 5.0f;
+                merged.push_back(r);
+            }
+    ProcessResult::SupercellClass sc;
+    sc.trend = MergedDoubledClassTrend(merged, *parity);
+    REQUIRE(sc.trend.size() >= 4);
+    CHECK(sc.trend.front().excess < 0.4);
+    CHECK(sc.trend.back().excess > 0.7);
+    CHECK(!DoubledCellRises(sc).empty());
 }
 
 TEST_CASE("FitDoubledClassTranslation", "[process]") {