diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 981e3b03f..a3882c0e7 100644 --- a/docs/RUGNUX_REPORT.md +++ b/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 (`_doubled.mtz`, or `_halfcell.mtz` where the doubled cell was adopted, with the rest of its files under the same diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index d660f7987..7ff390fca 100644 --- a/rugnux/ResultReport.cpp +++ b/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))); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 974069985..a9265b280 100644 --- a/rugnux/Rugnux.h +++ b/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 std::shared_ptr 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> DoubledClassParity(const UnitCell &doubled_cell, const UnitCell &half_cell, + char half_centering); +std::vector MergedDoubledClassTrend(const std::vector &merged, + const std::array &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 diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index d260bffab..cef22f28c 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1695,10 +1695,23 @@ std::string DoubledCellRises(const ProcessResult::SupercellClass &sc) { (outer->excess - low.excess) / se); } -DoubledClassTranslation FitDoubledClassTranslation(const std::vector &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 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> 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 parity{}; for (int c = 1; c < 8; c++) { const std::array 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 0.01) + return std::nullopt; + return parity; +} + +std::vector MergedDoubledClassTrend(const std::vector &merged, + const std::array &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 . 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 out; + std::array s_odd{}, s2_odd{}, s_even{}, s_ios{}; + std::array 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(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(n_odd[sh]), ne = static_cast(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 &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 {:.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 diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 7114e22f3..c6828a842 100644 --- a/tests/RugnuxTest.cpp +++ b/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{0, 0, 1})); + std::mt19937 rng(7); + std::exponential_distribution wilson(1.0); + std::vector 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(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(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]") {