From 8fa5b14bc03f97fa61f7d47b5cc7315cd412445f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 13 Oct 2026 19:42:25 +0200 Subject: [PATCH 1/3] rugnux: the rotation-scale walk stops where the next step would move the sweep end by under half a frame Each refit reads a share of the error left, so the walk's steps contract by about the ratio of the last two. Where the step that ratio predicts moves the end of the sweep by less than half a frame, the refit (an integration of the pre-pass frames and a post-refinement) and the indexing probe after it are not made; the scale the validation frames last preferred is adopted, as before. Measured on the walk sets of the battery against rc175 edde4439c: p.hkl and MTZ data identical; it fired on one open set (2 refits instead of 3, same adopted scale), about -7 s there. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- rugnux/Rugnux.cpp | 12 +++++++++++- rugnux/Rugnux.h | 3 ++- rugnux/RugnuxPasses.cpp | 3 ++- tests/RotationScaleWalkTest.cpp | 20 ++++++++++++++++++-- 4 files changed, 33 insertions(+), 5 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 179f7f311..9c4e6af05 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -134,7 +134,7 @@ bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e, double signi RotationScaleWalk WalkRotationScale(double first_fit, const std::function &index_at, const std::function(float)> &refit_at, - int max_rounds) { + int max_rounds, int64_t sweep_images) { RotationScaleWalk walk; auto next = static_cast(first_fit); if (next == walk.scale) @@ -145,13 +145,23 @@ RotationScaleWalk WalkRotationScale(double first_fit, }; walk.evidence = index_at(walk.scale); walk.trail = fmt::format("{:.5f}: {:.1f}%", walk.scale, score(walk.evidence)); + double last_step = 0.0; for (int round = 0; round < max_rounds && next != walk.scale; ++round) { const ValidationSpotEvidence e = index_at(next); walk.trail += fmt::format(", {:.5f}: {:.1f}%", next, score(e)); if (!ValidationEvidencePrefers(walk.evidence, e)) break; + const double step = next - walk.scale; walk.scale = next; walk.evidence = e; + // Each refit reads a share of the error left, so the steps contract by about the ratio of the + // last two. Where the step that ratio predicts moves the end of the sweep by less than half a + // frame, the scale it would reach indexes the frames as this one does, and the refit is not made. + if (last_step != 0.0 && std::fabs(step * step / last_step) * static_cast(sweep_images) < 0.5) { + walk.trail += fmt::format(" (next step predicted {:.5f}, under half a frame)", step * step / last_step); + break; + } + last_step = step; const auto fit = refit_at(walk.scale); if (!fit) break; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 5bba0695a..e1d703bba 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -284,6 +284,7 @@ std::pair ValidationEvidenceDifference(const ValidationSpotEvide // index_at(k): the validation evidence of the lattice indexed with the stored angles scaled by k // refit_at(k): the scale the post-refinement fits on reflections integrated at k, relative to k; // empty where it fitted none +// sweep_images: the frames of the sweep, against which a step too small to refit for is measured struct RotationScaleWalk { float scale = 1.0f; // adopted; 1 = the stored angles stand ValidationSpotEvidence evidence; // at the adopted scale @@ -292,7 +293,7 @@ struct RotationScaleWalk { RotationScaleWalk WalkRotationScale(double first_fit, const std::function &index_at, const std::function(float)> &refit_at, - int max_rounds); + int max_rounds, int64_t sweep_images); // The first pass found no lattice of the crystal in the frames - what a sample with no crystal in the // beam ends in. A statement about the data, not an error in the input: the CLI writes a report with diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 986c273f6..4369dd9af 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -729,7 +729,8 @@ RotationScaleWalk Rugnux::WalkRotationScaleBeside(const ProcessResult &pass1, Ro }; constexpr int MAX_ROTATION_SCALE_ROUNDS = 8; const RotationScaleWalk walk = WalkRotationScale(pass1.post_refine->rotation_scale, index_at, - refit_at, MAX_ROTATION_SCALE_ROUNDS); + refit_at, MAX_ROTATION_SCALE_ROUNDS, + experiment_.GetImageNum()); std::string text; if (!walk.trail.empty()) { text = fmt::format( diff --git a/tests/RotationScaleWalkTest.cpp b/tests/RotationScaleWalkTest.cpp index 617610485..7f942cca4 100644 --- a/tests/RotationScaleWalkTest.cpp +++ b/tests/RotationScaleWalkTest.cpp @@ -31,9 +31,9 @@ namespace { return 1.0 + 0.7 * (true_scale / k - 1.0); } - RotationScaleWalk Walk(double first_fit) { + RotationScaleWalk Walk(double first_fit, int64_t sweep_images = 3600) { return WalkRotationScale(first_fit, [this](float k) { return IndexAt(k); }, - [this](float k) { return RefitAt(k); }, 8); + [this](float k) { return RefitAt(k); }, 8, sweep_images); } }; } @@ -65,6 +65,22 @@ TEST_CASE("WalkRotationScale_ReachesTheFixedPoint", "[RotationScale]") { CHECK_FALSE(walk.trail.empty()); } +TEST_CASE("WalkRotationScale_StopsWhereTheNextStepIsUnderHalfAFrame", "[RotationScale]") { + // Spots enough for the validation frames to prefer every step: what ends the walk is the step the + // contraction predicts, measured against the sweep. Each refit reads 70 % of the error left, so the + // steps contract by 0.3. + SyntheticSweep short_sweep{0.97}; + short_sweep.spots = 1000000000; + const auto a = short_sweep.Walk(short_sweep.RefitAt(1.0f).value(), 360); + SyntheticSweep long_sweep{0.97}; + long_sweep.spots = 1000000000; + const auto b = long_sweep.Walk(long_sweep.RefitAt(1.0f).value(), 36000); + CHECK(short_sweep.refit_calls < long_sweep.refit_calls); + // Where it stops, what is left of the error moves the end of the sweep by under a frame. + CHECK(std::fabs(a.scale - 0.97) * 360 < 1.0); + CHECK(std::fabs(b.scale - 0.97) * 36000 < 1.0); +} + TEST_CASE("WalkRotationScale_StoredAnglesStand", "[RotationScale]") { SECTION("A healthy stage: a fit off by noise scores no better than the stored angles") { SyntheticSweep sweep{1.0}; From 7bef8d743c504060873018d4dbd27ac0a53f5e20 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 13 Oct 2026 19:42:25 +0200 Subject: [PATCH 2/3] rugnux: a canonical pass the geometry walk replaces stops after its space-group search arms The walk's indexing probes start as the canonical pass's post-refinement fit ends and run beside its merge and space-group search. The pass asked whether the walk takes a round only before the search and before its analyses; it now also asks after the search's two arms, before the space-group decisions and the merge in the adopted group, without waiting for anything. A pass the walk is going to replace stops there instead of finishing the search. (Waiting for the probes before the search was measured and dropped: it cost 0.2-3.6 s on sets where the walk took no round.) Exact by construction (only a pass that is thrown away is cut short); p.hkl and MTZ identical on the ci tier and 49 extra sets. Fired on one geometry-walk set, alternated runs 33.8 -> 31.6 s. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- rugnux/RugnuxScaleMerge.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 3e49e3c99..f609893f1 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1523,7 +1523,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (guard_beside) { guard_merge = guard_ahead.get(); guard_log.ReplayInto(logger); - record_search_merge(*guard_merge); + record_search_merge(*guard_merge, result); } if (!writer && GeometryWalkTakesRound(result, p.postrefine_ahead)) { logger.Info("The geometry walk will re-integrate at the geometry the post-refinement asks for, " @@ -1789,6 +1789,15 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { all_merge.reset(); rsm_all.reset(); + // The walk's indexing probes, which the fit starts as it ends, run beside the search; where + // their answer came in during its arms, a pass the walk replaces stops here, before the + // decisions and the merge in the adopted group. Nothing waits for it. + if (!writer && GeometryWalkTakesRound(result, p.postrefine_ahead)) { + logger.Info("The geometry walk will re-integrate at the geometry the post-refinement asks for, " + "so this pass stops after its space-group search arms"); + return true; + } + // Everything from here to the adopted group's merge: the lattice re-asks and the choices // between the candidates the arms left. TimingMark decisions_mark(logger, "space-group decisions"); From 94a51ea1ae88646f01eb5af88589abb9adc0561a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 13 Oct 2026 19:42:25 +0200 Subject: [PATCH 3/3] rugnux: the beam-centre arm at the measured centre integrates a block-strided third of the pre-pass first The measured-centre arm is a whole geometry pre-pass. It now integrates first a third of the pre-pass frames, in 3k + 1 blocks of about 10 deg over the whole pre-pass arc with every third kept, and the file's centre merges the same frames (its outcomes for the other frames are moved out while that engine ingests). Where the measured centre loses on that subset on every line MeasuredCentreWins draws by more than NEAR_TIE_SIGMA times the line's noise (Poisson for the strong counts, the Fisher standard error for a CC1/2), it has lost and no full arm is run. Anything else, a near tie or a win, puts the run back as the file's centre left it and runs the arm over all the pre-pass frames, deciding as before. The run made again (decide_all_late_) runs the arm on all frames. The same subset was tried for the metric-symmetry demoted arm and the walk refits and is not used there: the held-out residual's standard error does not grow on a third of the frames, so a subset margin overstates the evidence (one set's metric call flipped at 2.07 sigma), and a subset refit moved one walk to a worse fixed point. Against rc175 edde4439c, ci + 49 extra + 9 private sets: p.hkl and MTZ data identical on all 86. The arm lost on the subset on 22 sets (no full arm), was a tie/win on 9 (full arm as before, +0.3-1.7 s). Alternated runs: one set 62.1 -> 55.2 s, another 46.2 -> 42.4 s. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- rugnux/Rugnux.cpp | 16 ++++ rugnux/Rugnux.h | 37 +++++++++ rugnux/RugnuxImagePass.cpp | 5 ++ rugnux/RugnuxPasses.cpp | 155 ++++++++++++++++++++++++++---------- rugnux/RugnuxScaleMerge.cpp | 62 ++++++++++++--- tests/RugnuxTest.cpp | 30 +++++++ 6 files changed, 252 insertions(+), 53 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 9c4e6af05..7fd983129 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -177,6 +177,10 @@ namespace { // The least rotation the two-pass geometry pre-pass integrates, whatever --prepass-fraction asks. constexpr double PREPASS_MIN_ROTATION_DEG = 120.0; + + // About how wide a block of the pre-pass subset is (PrepassSubset). Wide against a reflection's + // rocking width, so that few of the reflections the subset holds are cut by a block's edges. + constexpr double PREPASS_SUBSET_BLOCK_DEG = 10.0; } namespace rugnux_internal { @@ -667,6 +671,18 @@ int Rugnux::PrepassEnd(int start_image, int end_image) const { return std::min(end, end_image); } +std::vector Rugnux::PrepassSubset(int prepass_images) const { + double arc_deg = 0.0; + if (const auto g = reader_.GetDataset()->experiment.GetGoniometer()) + arc_deg = std::abs(g->GetIncrement_deg()) * config_.stride * prepass_images; + const int k = std::max(1, static_cast(std::lround(arc_deg / (3.0 * PREPASS_SUBSET_BLOCK_DEG)))); + const int64_t blocks = 3 * k + 1; + std::vector keep(prepass_images); + for (int64_t o = 0; o < prepass_images; ++o) + keep[o] = (o * blocks / prepass_images) % 3 == 0; + return keep; +} + ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) { Logger logger("Rugnux"); ProcessResult result; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index e1d703bba..e4ed4d694 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -305,6 +305,16 @@ public: : JFJochException(JFJochExceptionCategory::InputParameterInvalid, r), reason(r) {} }; +// The numbers of a search merge (ProcessResult::search_merge_*) the beam-centre arms are judged on +// (MeasuredCentreWins), with the counts their noise is read from. +struct SearchMergeSignal { + int64_t strong = 0; // reflections at I/sigma >= 2 + int64_t reflections = 0; + double cc_half = NAN; + double inner_cc_half = NAN; // in the lowest-resolution shell ... + int64_t inner_reflections = 0; // ... of this many reflections +}; + struct ProcessResult { bool cancelled = false; uint64_t images_processed = 0; @@ -447,6 +457,12 @@ struct ProcessResult { // reaches. An empty reflection is as likely to land in one tail as in the other, so this is the // share of the strong count that is noise rather than signal (see RefinedPassIsWorse). int64_t search_merge_negative_reflections = 0; + // The reflections of that merge's lowest-resolution shell, which search_merge_inner_cc_half is over. + int64_t search_merge_inner_reflections = 0; + // The same search merge made of the pre-pass subset alone (Rugnux::PrepassSubset), on a geometry + // pre-pass followed by a beam-centre arm: that arm integrates those frames first, and the two + // centres are compared on them (see RunAllPasses). + std::optional prepass_subset_merge; // Per-reflection (I, sigma) of the final merged reflections, for the ISa diagnostic: I/sigma // plotted against I flattens off at the asymptote the error model reports as ISa, so the plot shows @@ -741,6 +757,14 @@ bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b); // geometry is the more accurate, answered on the signal each measured; where they differ it is which // metric the crystal has, answered on the search merges. std::string MeasuredCentreWins(const ProcessResult &file_arm, const ProcessResult &measured_arm); +// The search-merge numbers of a pass, as MeasuredCentreWins reads them. +SearchMergeSignal SearchSignal(const ProcessResult &r); +// Whether the measured centre loses on every line MeasuredCentreWins draws by more than `sigmas` +// times that line's noise - the Poisson noise of a count, the standard error of a CC1/2 - with both +// arms merged from the same frames. The run asks it at the near-tie distance of its first passes +// (FirstPassRun::NEAR_TIE_SIGMA); anything else, a near tie or a win, is decided on all the frames. +bool MeasuredCentreClearlyLoses(const SearchMergeSignal &file_arm, const SearchMergeSignal &measured_arm, + bool same_lattice, double sigmas); // What a cell refined UNDER a group's own constraint leaves of that group's metric violation: float // rounding on the lattice vectors, and nothing else (the constraint holds the angles at 90 exactly, @@ -890,6 +914,13 @@ class Rugnux { // Where the carried centre was found by reading less of each frame, the depth it was found at, // for the measured-centre arm to start from. std::optional beam_center_alternative_state_; + // Whether this geometry pre-pass is the beam-centre arm at the measured centre (RunAllPasses). The + // pre-pass at the file's centre merges the pre-pass subset as well where that arm follows it + // (ProcessResult::prepass_subset_merge); the arm itself does not. + bool measured_centre_arm_ = false; + // Whether a geometry pre-pass integrates only the pre-pass subset (PrepassSubset): the beam-centre + // arm at the measured centre, first. + bool prepass_subset_only_ = false; // The same pre-scan measurement, kept for the whole run instead of being consumed by the first // pass. Post-refinement measures how far it may move the beam from the nearer of the nominal @@ -1140,6 +1171,12 @@ class Rugnux { // sweep is still being collected, and an offline run makes them from the same frames. It is the // first config_.prepass_fraction of the sweep, but at least PREPASS_MIN_ROTATION_DEG (Rugnux.cpp). [[nodiscard]] int PrepassEnd(int start_image, int end_image) const; + // A third of the pre-pass frames in blocks spread over its whole arc: the arc is cut into 3k + 1 + // blocks of about PREPASS_SUBSET_BLOCK_DEG and every third is kept, the first and the last among + // them. What the beam-centre arm at the measured centre integrates first (prepass_subset_only_): + // its first pass and the angular reach of its post-refinement are a whole pre-pass's, and the two + // centres are merged from the same frames. Indexed by the pre-pass ordinal. + [[nodiscard]] std::vector PrepassSubset(int prepass_images) const; // The pre-scan's defective-pixel mask (HotPixelFinder), read on `sample` with the rings drawn // about `ring_centre` where the background measured one. frame(i, scratch) is sample frame i, as // the pre-scan kept it or read into scratch; null where it cannot be read. diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index 3a29de032..90060a378 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -228,6 +228,9 @@ void Rugnux::ImagePass(PipelineLocals &p) { } }; + // A pre-pass of the pre-pass subset (prepass_subset_only_) integrates its frames and no other. + const std::vector subset = prepass_subset_only_ ? PrepassSubset(prepass_images) : std::vector(); + auto full_worker = [&]() { pin_gpu(); // round-robin per worker thread; must precede engine construction // Before the analysis engine, which page-locks these bytes for its uploads and unregisters @@ -248,6 +251,8 @@ void Rugnux::ImagePass(PipelineLocals &p) { const int ordinal = next_ordinal.fetch_add(1); const int image_idx = start_image + ordinal * config_.stride; if (image_idx >= loop_end) break; + if (!subset.empty() && !subset[ordinal]) + continue; // Reuse the raw image frame to frame - a fresh 4-byte 16 Mpx frame is 72 MB of // untouched pages the decoder faults in one by one - except while the writer thread diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 4369dd9af..b620c9379 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -320,6 +320,32 @@ std::string MeasuredCentreWins(const ProcessResult &file_arm, const ProcessResul return {}; } +SearchMergeSignal SearchSignal(const ProcessResult &r) { + return {r.search_merge_strong_reflections, r.search_merge_reflections, r.search_merge_cc_half, + r.search_merge_inner_cc_half, r.search_merge_inner_reflections}; +} + +bool MeasuredCentreClearlyLoses(const SearchMergeSignal &file_arm, const SearchMergeSignal &measured_arm, + bool same_lattice, double sigmas) { + // A count is Poisson, as the first pass's near ties read counts of reflections; a CC1/2 over n + // reflections has the standard error of Fisher's transform, (1 - cc^2) / sqrt(n - 3). + const auto cc_se = [](double cc, int64_t n) { + return n > 3 ? (1.0 - cc * cc) / std::sqrt(static_cast(n - 3)) : NAN; + }; + const auto loses = [sigmas](double margin, double noise) { + return std::isfinite(margin) && std::isfinite(noise) && margin < -sigmas * noise; + }; + if (same_lattice) + return loses(static_cast(measured_arm.strong - file_arm.strong), + std::sqrt(static_cast(measured_arm.strong + file_arm.strong))) + && loses(measured_arm.inner_cc_half - file_arm.inner_cc_half - MAX_CC_HALF_LOSS, + std::hypot(cc_se(measured_arm.inner_cc_half, measured_arm.inner_reflections), + cc_se(file_arm.inner_cc_half, file_arm.inner_reflections))); + return loses(measured_arm.cc_half - file_arm.cc_half - MAX_CC_HALF_LOSS, + std::hypot(cc_se(measured_arm.cc_half, measured_arm.reflections), + cc_se(file_arm.cc_half, file_arm.reflections))); +} + ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { // Only the refined passes of the rotation two-pass are judged against a header-geometry pass, and // only while this holds one; a second Run() on the same object must not inherit the last one's. @@ -434,31 +460,81 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { const auto file_last_rotation_result = last_rotation_result_; const bool file_supercell_settled = supercell_settled_; - experiment_ = experiment_before_first_pass; - // The measured spot width is a property of the crystal and was measured once, before this - // arm; put it back on the experiment so this pass narrows to the pre-pass radius from the - // same place the first arm did and leaves the same widened radii behind for pass 2. - if (bragg_adaptive_) - experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_); - experiment_.BeamX_pxl(alt_center[0]).BeamY_pxl(alt_center[1]); - prepass_mosaicity_.clear(); - prepass_detector_geometry_.reset(); - prepass_result_.reset(); - + const auto file_alternative_state = beam_center_alternative_state_; + // Everything the arm leaves on the run put back as the file's centre left it. + const auto keep_file_centre = [&] { + experiment_ = file_experiment; + config_.spot_finding = file_spot_finding; + prepass_mosaicity_ = file_mosaicity; + prepass_detector_geometry_ = file_geometry; + prepass_result_ = file_result; + bragg_adaptive_ = file_bragg_adaptive; + bragg_before_adaptive_ = file_bragg_before_adaptive; + bkg_starved_fraction_ = file_bkg_starved_fraction; + widened_crowding_ = file_widened_crowding; + last_rotation_result_ = file_last_rotation_result; + supercell_settled_ = file_supercell_settled; + }; + const auto run_arm = [&](bool subset) { + experiment_ = experiment_before_first_pass; + // The measured spot width is a property of the crystal and was measured once, before + // this arm; put it back on the experiment so this pass narrows to the pre-pass radius + // from the same place the first arm did and leaves the same widened radii behind for + // pass 2. + if (bragg_adaptive_) + experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_); + experiment_.BeamX_pxl(alt_center[0]).BeamY_pxl(alt_center[1]); + prepass_mosaicity_.clear(); + prepass_detector_geometry_.reset(); + prepass_result_.reset(); + beam_center_alternative_state_ = file_alternative_state; + beam_center_alternative_.reset(); + beam_center_same_lattice_alternative_.reset(); + measured_centre_arm_ = true; + prepass_subset_only_ = subset; + ProcessResult r; + try { + r = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true); + } catch (const std::exception &e) { + // A hypothesis that does not survive its own pass is not a failure of the run: + // the file's centre has already produced a complete first pass, and this arm + // exists only to be compared against it. A machine that ran out of a resource + // is, though: the file's centre would meet it next. + measured_centre_arm_ = prepass_subset_only_ = false; + if (IsFatalResourceError(e)) throw; + logger.Warning("Beam centre check: the pass at the measured centre did not complete ({})", + e.what()); + } + measured_centre_arm_ = prepass_subset_only_ = false; + r.pass_number = 1; + r.pass_count = 2; + return r; + }; + // The arm integrates the pre-pass subset first, where the file's centre merged it as well + // (ProcessResult::prepass_subset_merge), and is judged there against the file's centre. An + // arm that loses there beyond the noise has lost; anything else is judged again with the + // arm integrated over all the pre-pass frames, as the file's centre was. ProcessResult alt; - try { - alt = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true); - } catch (const std::exception &e) { - // A hypothesis that does not survive its own pass is not a failure of the run: the - // file's centre has already produced a complete first pass, and this arm exists only - // to be compared against it. A machine that ran out of a resource is, though: the - // file's centre would meet it next. - if (IsFatalResourceError(e)) throw; - logger.Warning("Beam centre check: the pass at the measured centre did not complete " - "({}) - the run keeps the file's centre", e.what()); - } - alt.pass_number = 1; - alt.pass_count = 2; + bool lost_on_subset = false; + if (pass1.prepass_subset_merge) { + alt = run_arm(true); + const SearchMergeSignal &f = *pass1.prepass_subset_merge; + const SearchMergeSignal m = SearchSignal(alt); + lost_on_subset = !cancelled_ && alt.has_merge_statistics + && MeasuredCentreClearlyLoses(f, m, ArmsHoldSameLattice(pass1, alt), + FirstPassRun::NEAR_TIE_SIGMA); + logger.Info("Beam centre check: on a third of the pre-pass frames the measured centre merges {} " + "reflections at I/sigma >= 2 against {}, CC1/2 before corrections {:.3f} against " + "{:.3f}, {:.3f} against {:.3f} in the lowest-resolution shell - {}", m.strong, + f.strong, m.cc_half, f.cc_half, m.inner_cc_half, f.inner_cc_half, + lost_on_subset ? "it loses beyond the noise" + : "not a loss beyond the noise, so both are judged on all the frames"); + if (!lost_on_subset && !cancelled_) { + keep_file_centre(); + alt = run_arm(false); + } + } else + alt = run_arm(false); // The question is asked again of what the two arms ENDED on. The check opened it on the // lattices the standard-floor search returned, before the short-axis pass and the geometry // post-refinement; either can replace that lattice, and on a crystal with an axis under the @@ -480,7 +556,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { "I/sigma >= 2 before corrections at the measured centre against {} at the " "file's", pass1.rotation_lattice_type->centering, gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), - alt.search_merge_strong_reflections, pass1.search_merge_strong_reflections); + alt.search_merge_strong_reflections, + lost_on_subset ? pass1.prepass_subset_merge->strong : pass1.search_merge_strong_reflections); else if (alt.rotation_lattice_type && pass1.rotation_lattice_type) logger.Info("Beam centre check: after their whole first passes the arms still differ - " "{}-centred {} at the file's centre against {}-centred {} at the measured one", @@ -488,7 +565,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), alt.rotation_lattice_type->centering, gemmi::crystal_system_str(alt.rotation_lattice_type->crystal_system)); - const std::string measured_wins = cancelled_ ? std::string() : MeasuredCentreWins(pass1, alt); + const std::string measured_wins = cancelled_ || lost_on_subset ? std::string() + : MeasuredCentreWins(pass1, alt); if (!measured_wins.empty()) { const auto cell_text = [](const ProcessResult &r) { if (!r.consensus_cell) @@ -503,23 +581,16 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { note("beam centre: the measured centre adopted, it merges better"); pass1 = std::move(alt); } else { - logger.Info("Beam centre check: the measured centre does not win ({} against " - "{} reflections at I/sigma >= 2, CC1/2 before corrections {:.3f} against " - "{:.3f} on the search merge) - the run keeps the file's centre", - alt.search_merge_strong_reflections, pass1.search_merge_strong_reflections, - alt.search_merge_cc_half, pass1.search_merge_cc_half); + if (lost_on_subset) + logger.Info("Beam centre check: the measured centre does not win - the run keeps the file's centre"); + else + logger.Info("Beam centre check: the measured centre does not win ({} against " + "{} reflections at I/sigma >= 2, CC1/2 before corrections {:.3f} against " + "{:.3f} on the search merge) - the run keeps the file's centre", + alt.search_merge_strong_reflections, pass1.search_merge_strong_reflections, + alt.search_merge_cc_half, pass1.search_merge_cc_half); note("beam centre: the file's centre kept, the measured one does not merge better"); - experiment_ = file_experiment; - config_.spot_finding = file_spot_finding; - prepass_mosaicity_ = file_mosaicity; - prepass_detector_geometry_ = file_geometry; - prepass_result_ = file_result; - bragg_adaptive_ = file_bragg_adaptive; - bragg_before_adaptive_ = file_bragg_before_adaptive; - bkg_starved_fraction_ = file_bkg_starved_fraction; - widened_crowding_ = file_widened_crowding; - last_rotation_result_ = file_last_rotation_result; - supercell_settled_ = file_supercell_settled; + keep_file_centre(); } // The arm's own pass asks the same question again at its own centre; the run has already // decided it, so nothing is to consume a second answer. diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index f609893f1..a2b5f44a4 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -754,6 +754,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const bool is_rotation = experiment_.IsRotationIndexing(); // rotation indexing -> rotation scaling/merge std::optional rsm; std::optional guard_merge; // what the quality guard reads, where not the search merge + std::optional subset_merge; // the pre-pass subset's, for a beam-centre arm (below) std::optional prepass_postrefine_obs; // The canonical pass's post-refinement gather, run beside the chain (see there): set once it has // read the reflections, which nothing may change before. @@ -1097,6 +1098,38 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (merge_first_part) logger.Info("Two-pass: the quality guard reads this pass's merge of the first {} images, " "the ones the pre-pass integrated", prepass_images); + // A geometry pre-pass that a beam-centre arm follows merges the pre-pass subset as well - the + // frames that arm integrates first - for the two centres to be compared on the same frames + // (RunAllPasses). The other frames are taken out of the outcomes while its engine ingests, and + // it is merged before this pass's own engine is on the card. The run made again compares + // the two on all the frames. + if (geometry_prepass && beam_center_alternative_ && !measured_centre_arm_ && !decide_all_late_) { + TimingMark mark(logger, "search merge (pre-pass subset)"); + auto &outcomes = indexer->GetIntegrationOutcome(); + const std::vector subset = PrepassSubset(prepass_images); + std::vector held(outcomes.size()); + const auto swap_others = [&] { + for (size_t o = 0; o < outcomes.size(); ++o) + if (o >= subset.size() || !subset[o]) + std::swap(held[o], outcomes[o]); + }; + RotationScaleMerge subset_rsm(experiment_, outcomes, result.consensus_cell, + static_cast(config_.scaling_iter), config_.nthreads, logger, ""); + subset_rsm.SetWriteBackPerFrameScale(false); + swap_others(); + try { + subset_rsm.Ingest(); + } catch (...) { + swap_others(); + throw; + } + swap_others(); + subset_rsm.SetExportScaledFulls(false); + subset_rsm.SetFrenchWilson(false); + auto r = subset_rsm.Run(!experiment_.GetGemmiSpaceGroup().has_value(), + /*full_stats=*/false, /*measure_cc_before_corrections=*/true); + subset_merge = ScaleMergeResult{std::move(r.merged), std::move(r.statistics), r.cc_half_before_corrections}; + } // A reference MTZ is allowed for rotation: it fixes the space group / cell (on the CLI) and // resolves the indexing ambiguity (below), but is NOT used to scale - the rotation merge stays // self-consistent. @@ -1317,14 +1350,16 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // surfaces - where their final merges need not even be in the same group. See RunAllPasses. Where // the pre-pass integrated only part of the sweep, it is the same merge of those frames (above); // a guard merging beside the chain is read where the pass first needs it (below). - const auto record_search_merge = [&](const ScaleMergeResult &counted) { + const auto record_search_merge = [&](const ScaleMergeResult &counted, ProcessResult &into) { const auto &o = counted.statistics.overall; - result.search_merge_completeness_measured = o.possible_unique_reflections > 0; - result.search_merge_completeness = result.search_merge_completeness_measured + into.search_merge_completeness_measured = o.possible_unique_reflections > 0; + into.search_merge_completeness = into.search_merge_completeness_measured ? 100.0 * o.unique_reflections / o.possible_unique_reflections : 0.0; - result.search_merge_cc_half = counted.cc_half_before_corrections; - if (!counted.statistics.shells.empty()) - result.search_merge_inner_cc_half = counted.statistics.shells.front().cc_half; + into.search_merge_cc_half = counted.cc_half_before_corrections; + if (!counted.statistics.shells.empty()) { + into.search_merge_inner_cc_half = counted.statistics.shells.front().cc_half; + into.search_merge_inner_reflections = counted.statistics.shells.front().unique_reflections; + } // ... and the reflections it holds on the principal axial rows, which the guard needs // separately because CC1/2 cannot express them (see RunAllPasses). Only the low-order // part is counted: a screw is read off the strong start of its row - the evidence gate @@ -1352,7 +1387,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } }; constexpr float AXIAL_ROW_LOW_ORDER_D_A = 4.0f; - result.search_merge_axial_reflections = std::count_if( + into.search_merge_axial_reflections = std::count_if( counted.merged.begin(), counted.merged.end(), [&](const MergedReflection &r) { return ((r.h == 0) + (r.k == 0) + (r.l == 0)) == 2 && r.d > AXIAL_ROW_LOW_ORDER_D_A && centring_allows(r); @@ -1360,19 +1395,24 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // The signal the pass measured, and the size of the merge it measured it in, which is // what the guard decides on. constexpr float STRONG_I_OVER_SIGMA = 2.0f; - result.search_merge_strong_reflections = std::count_if( + into.search_merge_strong_reflections = std::count_if( counted.merged.begin(), counted.merged.end(), [&](const MergedReflection &r) { return r.sigma > 0 && r.I / r.sigma >= STRONG_I_OVER_SIGMA && centring_allows(r); }); - result.search_merge_negative_reflections = std::count_if( + into.search_merge_negative_reflections = std::count_if( counted.merged.begin(), counted.merged.end(), [&](const MergedReflection &r) { return r.sigma > 0 && r.I / r.sigma <= -STRONG_I_OVER_SIGMA && centring_allows(r); }); - result.search_merge_reflections = std::count_if(counted.merged.begin(), counted.merged.end(), + into.search_merge_reflections = std::count_if(counted.merged.begin(), counted.merged.end(), centring_allows); }; if (!guard_beside) - record_search_merge(guard_merge ? *guard_merge : sm); + record_search_merge(guard_merge ? *guard_merge : sm, result); + if (subset_merge) { + ProcessResult subset_view; + record_search_merge(*subset_merge, subset_view); + result.prepass_subset_merge = SearchSignal(subset_view); + } // Rotation two-pass geometry pre-pass: this first scale/merge has now fitted a frame-order-smoothed // mosaicity and written it back onto the per-frame outcomes. Capture it for the second pass's Bragg diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 48a3c37ba..bb8367cb7 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -529,3 +529,33 @@ TEST_CASE("MeasuredCentreWins", "[process]") { CHECK(!MeasuredCentreWins(file, arm(monoclinic, mono, 0.95, 9000)).empty()); CHECK(MeasuredCentreWins(file, arm(monoclinic, mono, 0.88, 20000)).empty()); } + +TEST_CASE("MeasuredCentreClearlyLoses", "[process]") { + // Same lattice: the measured centre loses only where it is behind on the signal by more than the + // Poisson noise of the two counts AND its inner shell does not agree better by the margin. + const SearchMergeSignal file{.strong = 9000, .reflections = 40000, .cc_half = 0.9, + .inner_cc_half = 0.98, .inner_reflections = 2000}; + SearchMergeSignal m = file; + m.strong = 8500; // 500 behind, noise sqrt(17500) = 132 + CHECK(MeasuredCentreClearlyLoses(file, m, true, 2.0)); + m.strong = 8800; // 200 behind: within two sigmas + CHECK_FALSE(MeasuredCentreClearlyLoses(file, m, true, 2.0)); + // A file centre whose few inner reflections disagree: the measured centre's inner shell might + // still carry it, so behind on the signal is not yet a loss. + SearchMergeSignal inconsistent = file; + inconsistent.inner_cc_half = 0.6; + inconsistent.inner_reflections = 20; + m.strong = 8500; + m.inner_cc_half = 0.64; + m.inner_reflections = 20; + CHECK_FALSE(MeasuredCentreClearlyLoses(inconsistent, m, true, 2.0)); + + // Different lattices: the CC1/2 line, 0.05 above the file's centre. + SearchMergeSignal d = file; + d.cc_half = 0.80; + CHECK(MeasuredCentreClearlyLoses(file, d, false, 2.0)); + d.cc_half = 0.95; // exactly on the line + CHECK_FALSE(MeasuredCentreClearlyLoses(file, d, false, 2.0)); + // Nothing measured is no loss: it is decided on all the frames. + CHECK_FALSE(MeasuredCentreClearlyLoses(file, SearchMergeSignal{}, false, 2.0)); +}