diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 179f7f311..7fd983129 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; @@ -167,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 { @@ -657,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 5bba0695a..e4ed4d694 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 @@ -304,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; @@ -446,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 @@ -740,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, @@ -889,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 @@ -1139,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 986c273f6..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. @@ -729,7 +800,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/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 3e49e3c99..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 @@ -1523,7 +1563,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 +1829,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"); 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}; 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)); +}