diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index da2765a2a..7eb93b418 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/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(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(row_present_esq.at(a.row).size()); } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index e4201767b..21227a34b 100644 --- a/rugnux/Rugnux.h +++ b/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 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 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 parity{}; + std::array eps{}; + std::vector shells; +}; +DoubledClassTranslation FitDoubledClassTranslation(const std::vector &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 force_rotation_result_; // The indexing result the last full pass integrated on, which the doubled-cell hypothesis doubles. std::optional 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, diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index f3d7b8091..923acff63 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1415,27 +1415,24 @@ std::optional 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 &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 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(p[0]) + v[1] * static_cast(p[1]) + v[2] * static_cast(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 sum{}, sum2{}; + std::array 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(r.I) * r.I; cnt[sh]++; + } + } + std::array 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(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 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(r.h), static_cast(r.k), static_cast(r.l), + r.I / (sum[sh] / static_cast(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 &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 best{}; + double best_cc = -INFINITY; + const size_t coarse_stride = std::max(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 centre = best; + const int nstep = static_cast(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 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(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 _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 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(); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 4b3b058c0..ac30c40e7 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/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(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) { diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 60f00785b..7e7dbdaff 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include #include @@ -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 wilson(1.0); + std::vector 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(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(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{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) {