From 2e0a3bb2553b2aa4787b00ddb62cd6ccecea8177 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 7 Oct 2026 15:48:43 +0200 Subject: [PATCH] rugnux: move RunAllPasses and the pass arbitration to RugnuxPasses.cpp A verbatim move of RunAllPasses, RefinedPassIsWorse, ArmsHoldSameLattice and MeasuredCentreWins, with the helpers only they read (MAX_CC_HALF_LOSS, SearchMergeCompletenessText, PrimitiveCellVolume, SignalText, SameSetting) in the new file's anonymous namespace. The comment describing RefinedPassIsWorse, which sat at the end of that anonymous namespace, now sits on the function. TwinZoneVerdict stays behind for the scale/merge phase. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/CMakeLists.txt | 1 + rugnux/Rugnux.cpp | 1152 ------------------------------------ rugnux/RugnuxPasses.cpp | 1242 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 1243 insertions(+), 1152 deletions(-) create mode 100644 rugnux/RugnuxPasses.cpp diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index eb8595bba..91e575e14 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -7,6 +7,7 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxStills.cpp RugnuxCalibration.cpp RugnuxPreScan.cpp + RugnuxPasses.cpp RugnuxDefaults.cpp RugnuxDefaults.h Rugnux.h diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 73f7df09c..9fdfe60dc 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -512,88 +512,6 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { } namespace { - // How much better the alternative beam centre's search merge has to be for the run to move off - // the centre the file gave it (the beam-centre arbitration in RunAllPasses). - constexpr double MAX_CC_HALF_LOSS = 0.05; - - // Say so rather than printing a zero that reads as a measurement. ("not compared", not "not - // measured": completeness IS measured and reported; what did not happen is its use in arbitrating - // between the two passes.) - std::string SearchMergeCompletenessText(const ProcessResult &pass1, const ProcessResult &pass2) { - if (!pass1.search_merge_completeness_measured || !pass2.search_merge_completeness_measured) - return "completeness not compared"; - return fmt::format("completeness {:.1f}% vs {:.1f}%", pass2.search_merge_completeness, - pass1.search_merge_completeness); - } - - // The volume of the pass's primitive cell, 0 where it has no cell. - double PrimitiveCellVolume(const ProcessResult &r) { - if (!r.consensus_cell) - return 0.0; - return std::abs(CrystalLattice(*r.consensus_cell).ToPrimitive(r.consensus_centering.value_or('P')) - .CalcVolume()); - } - - // The signal the two passes measured, as the quality guard reads it (RefinedPassIsWorse). - std::string SignalText(const ProcessResult &pass1, const ProcessResult &pass2) { - return fmt::format("{} reflections at I/sigma >= 2 less the {} at <= -2 before corrections, " - "against {} less {}", pass2.search_merge_strong_reflections, - pass2.search_merge_negative_reflections, pass1.search_merge_strong_reflections, - pass1.search_merge_negative_reflections); - } - - // Whether the two passes describe their lattice in the same setting, which is what makes a count - // of reflections at a given Miller index the same count on both sides. A pass that re-indexes at - // the refined geometry can settle on a centred setting of the same lattice the other pass held - // primitive; the cells are then different bases and only the counts taken over all reflections - // survive the change. - bool SameSetting(const ProcessResult &pass1, const ProcessResult &pass2) { - if (pass1.consensus_centering != pass2.consensus_centering - || !pass1.consensus_cell || !pass2.consensus_cell) - return false; - const UnitCell &c1 = *pass1.consensus_cell, &c2 = *pass2.consensus_cell; - const auto same_length = [](double a, double b) { return std::abs(a - b) <= 0.05 * (a + b) / 2; }; - return same_length(c1.a, c2.a) && same_length(c1.b, c2.b) && same_length(c1.c, c2.c); - } - - // Pass 2 is normally the better answer, which is why it is the canonical output - but it is not - // guaranteed to be. Returns why the refined pass is worse, in the words the run reports it in, or - // an empty string when it is not. - // - // The comparison is made on each pass's SEARCH merge, not on the final one. Each pass determines - // its own space group, so the two final merges can be in different groups, and a completeness or a - // correlation in different groups are not the same measurement. The search merge is P1 on both - // sides (or the group the user fixed), over the full resolution range and before any correction - // surface is fitted - and it is counted over the centring-allowed reflections on both sides, so - // that a pass which finds a centring the other did not is not judged on the extinguished half it - // merged as noise (see the capture in RunPipeline). - // - // What the two passes are compared on there is the SIGNAL each measured - the reflections it - // merged at I/sigma >= 2 - and not, as this guard decided until it was measured, that merge's - // CC1/2. Two measurements settled it: - // - Against an EXTERNAL arbiter (the site's own XDS merge of the same images, scored over a - // resolution range common to the arms) the count of unique reflections at I/sigma >= 2 agreed - // with which geometry is the more accurate one on 23 of 27 arm-dataset pairs. R_meas, CC1/2 - // and ISa - which measure whether a merge agrees with ITSELF - agreed on 13 of 27, a coin flip. - // - This merge's CC1/2 is pooled over the whole range, uncut and uncorrected, so the shells with - // no signal in them dominate it and they are exactly the shells a geometry move disturbs: it - // reads 0.13 on a crystal whose data merge at 0.995, and 0.93 against 0.76 on two runs whose - // merged data agree to 0.7% of R_meas. Deciding on it sent two crystals back to a geometry - // that ships 22% fewer reflections than the one they refused. - // - // Three ways the refined pass can be worse, then: it merges more unique reflections than the cell - // it settled on can hold, it measured decisively less signal, or it lost the axial rows. - // - // The completeness arm only bites when the search merge is in a fixed group (-S). De novo that - // merge is in P1 and does not count its possible reflections at all, so there is no completeness - // and it does not arm. Counting them there was measured and not taken: it costs 5.2 ms, but a - // merge that fills its own asymmetric unit reads 27-43% of the P1 hemisphere over 28 rotation - // crystals, so the 100.5% bound is never approached and no decision changes - it would only put a - // figure into the report that reads as the dataset's completeness and is not one. - // - // Asked twice, from the same numbers: in RunAllPasses, where the answer is acted on, and inside a - // refined pass as soon as its own search merge is in - a pass that already knows it will be thrown - // away leaves its report and its files unwritten (see the superseded flag in RunPipeline). // The twin-immune zone of the operators `higher` adds over `lower`, read on a P1 merge, as a // verdict: +1 the added operators are real (the zone is centric with margin), -1 they are a twin // law or a pseudo-symmetry (acentric with margin), 0 undecided - fewer than two orders apart @@ -666,1076 +584,6 @@ namespace { } -std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &pass2) { - constexpr double MAX_CREDIBLE_COMPLETENESS = 100.5; // rounding headroom, nothing more - // How much of the signal the refined pass has to lose before the run moves back off the - // geometry it measured. This is not a quality bar: it is set above what the two passes differ - // by for reasons that are not the geometry, and they never integrate quite the same way - the - // canonical pass integrates at the measured spot width where the pre-pass integrated at the - // default, and it predicts with the pre-pass's smoothed mosaicity. Measured over ~100 rotation - // crystals, that alone moves the count by -2% to +35%, and on a crystal where the two passes - // ran at the SAME geometry it moved it +10% (and the CC1/2 this guard used to decide on by - // -0.25). Below this, the difference says nothing about the geometry, and a comparison that - // cannot tell keeps the geometry the run measured. - constexpr double DECISIVE_SIGNAL_LOSS = 0.10; - // The signal is the count of reflections at I/sigma >= 2 LESS the count at I/sigma <= -2. An - // empty reflection - one the lattice predicts and the crystal does not diffract into - is as - // likely to land in either tail, so the negative tail is the merge's own measurement of how - // much of its strong count is noise, and what is left is what the crystal put there. It is what - // makes the two counts comparable when the passes hold different lattices: a pass on an n-fold - // supercell merges n times the reflections, n - 1 of every n of them empty, and their noise alone - // reached the strong count - measured, a header pass on a 7x supercell out-counted the refined - // pass on the crystal's own lattice and the run went back to the supercell, while that merge - // held more reflections at <= -2 than at >= 2. On a real superstructure the extra classes hold - // signal, and it is counted. - const auto signal = [](const ProcessResult &r) { - return static_cast(r.search_merge_strong_reflections - r.search_merge_negative_reflections); - }; - // ... and the loss has to be a loss rather than a dilution. The wider disk pulls weaker - // reflections into the merge and the error model is then refitted over them, so a pass can hold - // MORE reflections and fewer of them at I/sigma >= 2 without having measured less of the - // crystal: measured, one crystal merges 2.1x the reflections in its refined pass, a third of - // the strong ones, and its refined pass is the better answer by every final statistic. Where - // the refined pass covers no more of reciprocal space and measured less signal, the loss is a - // loss. Between two lattices the raw counts are no comparison: a lattice of n times the primitive - // volume holds n times the reflections to the same resolution, so against a header pass on a - // supercell the refined pass always held fewer and the clause never protected it. The count of - // the smaller lattice is scaled by the index n, which is an integer - the same lattice, whose - // volume two refinements put a percent or two apart, is index 1 and compares as it always did. - const auto coverage_does_not_grow = [&] { - double n1 = static_cast(pass1.search_merge_reflections); - double n2 = static_cast(pass2.search_merge_reflections); - const double v1 = PrimitiveCellVolume(pass1), v2 = PrimitiveCellVolume(pass2); - if (v1 > 0.0 && v2 > 0.0) { - if (v1 >= v2) - n2 *= std::round(v1 / v2); - else - n1 *= std::round(v2 / v1); - } - return n2 <= n1; - }; - - // The other arm: the axial rows. A bulk statistic is an average over tens of thousands of - // reflections and the principal axial rows are a few dozen of them, so a pass can measure half - // the row and move it by nothing - yet those are the reflections the systematic absences are - // read off, and losing them costs a screw axis and with it the space group. It happens because - // a reflection on an axial row crosses the Ewald sphere over a narrow range of the sweep, where - // a general reflection is spread over it: a pass that loses a wedge of the sweep loses whole - // axial rows while its bulk statistics improve. Measured over a dozen rotation crystals the two - // passes agree on this count to within 9% and usually exactly; the crystal this arm was written - // for measured 44%. - // - // It is asked only where the two counts are counts of the SAME reflections: the a and b of a - // C-centred setting are the face diagonals of the primitive cell, so "h00" is a different - // direction in the two merges and the two rows are not comparable at all. And it needs pass 1 - // to have measured enough of the rows for the difference to mean anything - below eight there - // is no zone a screw could be claimed on either side. A pass that measured decisively more - // signal has said something about the bulk of the data that a few dozen reflections cannot - // answer, and is kept. - constexpr double MIN_AXIAL_ROWS_KEPT = 0.75; - constexpr int64_t MIN_AXIAL_ROWS_TO_COMPARE = 8; - const double signal1 = signal(pass1); - const double signal2 = signal(pass2); - const int64_t axial1 = pass1.search_merge_axial_reflections; - const int64_t axial2 = pass2.search_merge_axial_reflections; - const bool measured_less = signal1 > 0 && signal2 < (1.0 - DECISIVE_SIGNAL_LOSS) * signal1 - && coverage_does_not_grow(); - const bool measured_more = signal1 > 0 && signal2 > (1.0 + DECISIVE_SIGNAL_LOSS) * signal1; - const bool lost_axial_rows = SameSetting(pass1, pass2) - && axial1 >= MIN_AXIAL_ROWS_TO_COMPARE - && axial2 < MIN_AXIAL_ROWS_KEPT * static_cast(axial1) - && !measured_more; - const bool worse = pass2.lattice_conflicts_with_fixed_sg - || (pass1.search_merge_completeness_measured - && pass2.search_merge_completeness_measured - && pass2.search_merge_completeness > MAX_CREDIBLE_COMPLETENESS) - || measured_less || lost_axial_rows; - if (!worse) - return {}; - return lost_axial_rows && !measured_less - ? fmt::format("it holds {} of the {} low-order axial reflections the header-geometry " - "pass measured, and its {} do not say it is better", axial2, axial1, - SignalText(pass1, pass2)) - : fmt::format("{}, {}", SearchMergeCompletenessText(pass1, pass2), SignalText(pass1, pass2)); -} - -bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b) { - if (!a.rotation_lattice_type || !b.rotation_lattice_type || !a.consensus_cell || !b.consensus_cell) - return false; - const auto primitive_volume = [](const ProcessResult &r) { - return std::abs(CrystalLattice(*r.consensus_cell) - .ToPrimitive(r.consensus_centering.value_or('P')).CalcVolume()); - }; - const double va = primitive_volume(a), vb = primitive_volume(b); - return a.rotation_lattice_type->crystal_system == b.rotation_lattice_type->crystal_system - && a.rotation_lattice_type->centering == b.rotation_lattice_type->centering - && std::max(va, vb) < 1.02 * std::min(va, vb); -} - -std::string MeasuredCentreWins(const ProcessResult &file_arm, const ProcessResult &measured_arm) { - if (!file_arm.has_merge_statistics || !measured_arm.has_merge_statistics) - return {}; - // Same lattice: which geometry is the more accurate, read off the signal each arm measured - the - // statistic the two-pass quality guard found to agree with an external arbiter where CC1/2 did not - // (RefinedPassIsWorse). Without its margin: that margin covers what the header and refined passes - // differ by for reasons that are not the geometry - a different spot width and mosaicity - and the - // two arms here integrate with the same settings, so the geometry is the only difference, and the - // file's centre has no standing the measured one lacks - neither indexed the whole spot list. - if (ArmsHoldSameLattice(file_arm, measured_arm)) { - if (measured_arm.search_merge_strong_reflections > file_arm.search_merge_strong_reflections) - return fmt::format("the measured centre measures more of the crystal - {} reflections at " - "I/sigma >= 2 before corrections against {}", - measured_arm.search_merge_strong_reflections, - file_arm.search_merge_strong_reflections); - // ...or its strongest reflections agree with themselves where the file's do not. The pooled - // CC1/2 is no arbiter here - the shells with no signal in them dominate it - but in the - // lowest-resolution shell every reflection is strong, and half-sets that disagree there are - // equivalents given different indices: a centre about one reflection spacing out. Measured, - // such a file centre merged 2 % MORE reflections at I/sigma >= 2 than the right one, at an - // inner-shell CC1/2 of 0.610 against 0.997. - if (measured_arm.search_merge_inner_cc_half > file_arm.search_merge_inner_cc_half + MAX_CC_HALF_LOSS) - return fmt::format("the measured centre's strongest reflections agree where the file's do not - " - "CC1/2 {:.3f} against {:.3f} in the lowest-resolution shell of the search merge", - measured_arm.search_merge_inner_cc_half, file_arm.search_merge_inner_cc_half); - return {}; - } - if (measured_arm.search_merge_cc_half > file_arm.search_merge_cc_half + MAX_CC_HALF_LOSS) - return fmt::format("the measured centre merges better - CC1/2 before corrections {:.3f} against " - "{:.3f} on the search merge", measured_arm.search_merge_cc_half, - file_arm.search_merge_cc_half); - return {}; -} - -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. - quality_guard_pass1_.reset(); - // Rotation two-pass geometry post-refinement: the first pass integrates at the header geometry and - // post-refines the detector geometry (distance + beam from the observed spot positions, cell scale + axis - // from phi_obs); the second pass re-indexes and re-integrates with the refined geometry. Only the second - // (refined-geometry) pass writes merged files, under the plain "_*" name. The first pass still runs - // under a "_01" prefix so that the _process.h5 it writes, when one is asked for, does not collide - // with the second pass's - the ice-ring flags and geometry stored in each are the ones that pass used. - if (config_.mode == ProcessMode::FullAnalysis && config_.rotation_postrefine_geometry - && experiment_.IsRotationIndexing()) { - Logger logger("Rugnux"); - // The checks run between the passes are steps of their own on the console; their probe passes - // post none. - const auto step = [observer](const std::string &title) { if (observer) observer->OnStep(title); }; - const auto note = [observer](const std::string &text) { if (observer) observer->OnNote(text); }; - const std::string base_prefix = config_.output_prefix; - const auto gonio_snapshot = experiment_.GetGoniometer(); - // What the first pass starts from, so a second first pass at a different beam centre can start - // from the same place (the beam-centre arbitration below). - const DiffractionExperiment experiment_before_first_pass = experiment_; - const int64_t max_spot_count_snapshot = experiment_.GetDatasetSettings().GetMaxSpotCount(); - header_distance_mm_ = experiment_.GetDetectorDistance_mm(); - prepass_detector_geometry_.reset(); - prepass_rotation_scale_.reset(); - prepass_result_.reset(); - force_rotation_result_.reset(); - bragg_adaptive_.reset(); - beam_center_same_lattice_alternative_.reset(); - - logger.Info("Rotation two-pass geometry post-refinement: first pass (header geometry) -> {}_01_*", - base_prefix); - // Keep an empty prefix empty (the "compute stats, persist nothing" mode): the "_01" suffix would - // make it non-empty and RunPipeline would then write stray _01_* files despite no output wanted. - config_.output_prefix = base_prefix.empty() ? base_prefix : base_prefix + "_01"; - auto pass1 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true); - pass1.pass_number = 1; - pass1.pass_count = 2; - if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } - - // Beam centre as a LATTICE hypothesis, decided on the merge. The first pass indexed both - // centres, found the same cell with a different metric symmetry read off it, and measured that - // the centre disagreement is by itself large enough to have decided that symmetry (see the - // beam-centre check in RunPipeline). The class is handed to GeometryRefiner, which holds the - // cell to it for the rest of the run, so the wrong one costs an axis - refined, predicted, - // integrated and merged with an equality the crystal does not have. - // - // The disagreement is two-sided, so there is no centre to prefer: over the corpus runs that - // reach here the file's centre finds the higher class about as often as the measured one does. - // Nothing before integration separates them either - both index the same frames, and the frame - // count is the one statistic the check already refuses to arbitrate on. So run the first pass - // again at the other centre and keep whichever arm MERGES better, each judged on its own search - // merge: P1, whole range, before any correction surface, which is the merge the two passes are - // already compared in at the end of this function. Because the alternative arm is a whole first - // pass it gets the short-axis pass and the geometry post-refinement of its own, which the - // question needs - on a crystal with an axis under the FFT floor the check runs before that - // pass and compares supercells. - // - // Adopting a different metric moves the run off the geometry it was given, so it has to be - // decisively better rather than merely ahead. - // - // The same two arms also answer a second question the first pass cannot: where neither centre - // indexes the whole spot list and both index the same lattice once the pass reads less of each - // frame, which centre is right (see the beam-centre check). The measured arm then starts from - // the depth the check found it at. - // - // And a third: where both centres index the same lattice in the same class but lie further - // apart than the geometry absorbs, the centre can be off by about one reflection spacing - the - // same cell, every reflection one index out over part of the sweep - and only the merge shows - // it. Running the first pass twice on every such run would cost a pass on half of them, so it - // is run where the file's centre merged inconsistently: the lowest-resolution shell of the - // search merge, where every reflection is strong and noise cannot pull the correlation down, - // reads its half-sets as disagreeing. Measured on a pseudo-tetragonal crystal whose two centres - // lie 10.5 px apart, one c* spacing being 11.3 px: 0.610 at the file's centre and 0.997 at the - // measured one, R_meas of the final merge 0.58 against 0.13. - constexpr double MIN_INNER_CC_HALF = 0.9; - if (!cancelled_ && !beam_center_alternative_ && beam_center_same_lattice_alternative_ - && pass1.search_merge_inner_cc_half < MIN_INNER_CC_HALF) { - logger.Warning("Beam centre check: the file's centre merges inconsistently - CC1/2 {:.3f} in " - "the lowest-resolution shell of the search merge, where every reflection is " - "strong - and the measured centre, which indexes the same lattice, is further " - "away than this geometry absorbs, so the run tries it", - pass1.search_merge_inner_cc_half); - beam_center_alternative_ = beam_center_same_lattice_alternative_; - } else if (beam_center_same_lattice_alternative_) - logger.Info("Beam centre check: the file's centre merges consistently - CC1/2 {:.3f} in the " - "lowest-resolution shell of the search merge - so the measured centre, which " - "indexes the same lattice, is not tried", pass1.search_merge_inner_cc_half); - beam_center_same_lattice_alternative_.reset(); - if (!cancelled_ && beam_center_alternative_) { - const auto alt_center = *beam_center_alternative_; - beam_center_alternative_.reset(); - logger.Info("Beam centre check: running the first pass again at the measured centre " - "({:.2f},{:.2f}) and judging the two on their merges", alt_center[0], - alt_center[1]); - note(fmt::format("the first pass runs again at the measured beam centre ({:.1f}, {:.1f})", - alt_center[0], alt_center[1])); - const DiffractionExperiment file_experiment = experiment_; - const auto file_spot_finding = config_.spot_finding; - const auto file_mosaicity = prepass_mosaicity_; - const auto file_geometry = prepass_detector_geometry_; - const auto file_result = prepass_result_; - - 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(); - - 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. - 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; - // 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 - // FFT floor the check compared two supercells the short-axis pass then threw away. Where - // both arms now hold the same lattice in the same class, the centre did not decide the - // symmetry after all: the two merges are two integrations of ONE hypothesis, their CC1/2 - // difference is the noise between them, and it is not asked. Same lattice means what it - // means in the check: same class and primitive volumes within 2 %. - // - // What is left is which of the two GEOMETRIES is the more accurate, answered on the signal - // each arm measured (MeasuredCentreWins). The post-refinement cannot be left to settle it: - // it moves the beam by what the spots ask of a lattice already fitted at that centre, which - // on the crystal that showed it was 0.15 px of a 12.6 px error. - const bool arms_agree = ArmsHoldSameLattice(pass1, alt); - if (arms_agree) - logger.Info("Beam centre check: after their whole first passes both centres hold the " - "same lattice - {}-centred {} - so the centre did not decide the metric " - "symmetry; the two are judged on the signal each measured, {} reflections at " - "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); - 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", - pass1.rotation_lattice_type->centering, - 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); - if (!measured_wins.empty()) { - const auto cell_text = [](const ProcessResult &r) { - if (!r.consensus_cell) - return std::string("no cell"); - const auto &c = *r.consensus_cell; - return fmt::format("{:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}", - c.a, c.b, c.c, c.alpha, c.beta, c.gamma); - }; - logger.Warning("Beam centre check: {}, on {} against {} - so the run adopts the measured " - "centre and the lattice it finds", measured_wins, cell_text(alt), - cell_text(pass1)); - 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); - 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; - } - // 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. - beam_center_alternative_.reset(); - beam_center_alternative_state_.reset(); - beam_center_same_lattice_alternative_.reset(); - if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } - } - - // Undo the pre-pass goniometer shift (a sub-range or strided run moves start/increment inside - // the pass), but KEEP an axis direction the pass changed. The rotation-axis sign rescue runs - // inside the pass and adopts the opposite axis when the file's own indexes almost nothing; - // restoring the snapshot wholesale threw that away, handing pass 2 the sign that failed and - // then applying pass 1's fitted rotation scale to the axis it was not fitted on. It - // self-healed - pass 2's rescue fires again and a wrong sign indexes almost no frames - but - // only by repeating a whole first pass, and only while "a wrong sign always scores below - // half" holds, which is an invariant nothing states or tests. - if (gonio_snapshot) { - GoniometerAxis restored = *gonio_snapshot; - if (const auto after_pass1 = experiment_.GetGoniometer()) - restored.Axis(after_pass1->GetAxis()); - experiment_.Goniometer(restored); - } - // A measured spot budget is the pass's own, not the run's: give the second pass the same list to - // measure from, so the two passes cannot ratchet each other down. The canonical pass then reads - // the budget off the refined geometry, where the spots that do lie on the lattice actually do. - experiment_.MaxSpotCount(max_spot_count_snapshot); - - // Apply the post-refined detector geometry for the second pass, keeping the header geometry so - // the run can go back to it if the refined pass turns out worse (see the quality guard below). - const std::array header_geometry{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), - experiment_.GetDetectorDistance_mm(), - experiment_.GetPoniRot1_rad(), - experiment_.GetPoniRot2_rad()}; - if (prepass_detector_geometry_) { - const auto &g = *prepass_detector_geometry_; - experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) - .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); - } - // The probe passes below and the lattice arms after them exist ONLY to measure - see there - so - // the count of passes the run made has to include them even though they write nothing. - int arm_passes = 0; - - // The goniometer rotation scale: the angles in the file are the commanded ones, so a stage that - // turned at the wrong rate is a geometry error like any other - but one pass 1's fit reads only - // in part (see WalkRotationScale), so it is walked to the fit's fixed point here and adopted - // only where the validation frames of the whole sweep prefer it. The walk runs at the detector - // geometry the second pass will, so every hypothesis, the stored angles included, is scored - // there, and at that geometry only. Its passes - // only measure - an indexing probe stops once the lattice is scored, a refit stops once the - // post-refinement has measured - and they leave nothing behind: the experiment, the geometry - // the second pass is to run at and pass 1's mosaicity (a width in degrees fitted against the - // stored angles) are put back when the walk ends. - std::string rotation_scale_walk; - if (!cancelled_ && gonio_snapshot && pass1.post_refine && !config_.rotation_scale) { - step("Rotation-scale check"); - const DiffractionExperiment before_walk = experiment_; - const auto geometry_before_walk = prepass_detector_geometry_; - const auto mosaicity_before_walk = prepass_mosaicity_; - // A probe pass at scale k, run on r - this run, or a copy of it (see below). - const auto probe_on = [&](Rugnux &r, RugnuxObserver *obs, float k, bool index_only) { - r.experiment_ = before_walk; - r.experiment_.Goniometer(ScaleRotation(*before_walk.GetGoniometer(), k)); - r.prepass_rotation_scale_ = k; - r.prepass_mosaicity_.clear(); - r.indexing_probe_only_ = index_only; - r.postrefine_probe_ = r.postrefine_probe_only_ = !index_only; - ProcessResult res; - try { - res = r.RunPipeline(obs, /*write_output=*/false, /*geometry_prepass=*/false); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - // Angles under which nothing indexes have scored nothing, which is what an empty - // result reads as. - logger.Info("Rotation scale {:.5f}: the probe pass did not complete ({})", k, e.what()); - } - r.indexing_probe_only_ = r.postrefine_probe_ = r.postrefine_probe_only_ = false; - return res; - }; - const auto probe = [&](float k, bool index_only) { - ++arm_passes; - return probe_on(*this, observer, k, index_only); - }; - // The walk always indexes at the stored angles first and at the pass-1 fit next (see - // WalkRotationScale). Both are known now and each is its own pass, so the one at the fit - // runs at the same time, on a copy of the run as it stands before either. What a probe - // leaves behind for the next pass - the spot-finding settings a first-pass ladder rung - // adopts; the experiment is put back by the probe itself - is checked when the fit's turn - // comes: if the probe at the stored angles changed it, the copy started from a state the - // sequential walk would not have, and the probe is run again here. - const float first_fit = static_cast(pass1.post_refine->rotation_scale); - std::shared_ptr fit_copy; - std::future at_fit; - const SpotFindingSettings walk_start_spot_finding = config_.spot_finding; - if (first_fit != 1.0f && !cancelled_) { - fit_copy = std::make_shared(*this); - at_fit = std::async(std::launch::async, [&probe_on, fit_copy, first_fit] { - return probe_on(*fit_copy, nullptr, first_fit, true); - }); - } - const auto index_at = [&](float k) { - if (at_fit.valid() && k == first_fit) { - ProcessResult r = at_fit.get(); - const bool same_start = - config_.spot_finding.high_resolution_limit == walk_start_spot_finding.high_resolution_limit - && config_.spot_finding.measured_ring_q_recipA == walk_start_spot_finding.measured_ring_q_recipA; - if (same_start) { - ++arm_passes; - config_.spot_finding = fit_copy->config_.spot_finding; - for (const auto &m : fit_copy->first_pass_memo_) - AddFirstPassMemo(m); - return r.validation_evidence; - } - } - return probe(k, true).validation_evidence; - }; - const auto refit_at = [&](float k) -> std::optional { - const ProcessResult r = probe(k, false); - if (!r.post_refine) - return std::nullopt; - logger.Info("Rotation scale {:.5f}: the post-refinement there fits {:.5f} on top of it, " - "held-out residual {:.3e}", k, r.post_refine->rotation_scale, - r.post_refine->held_out_before); - return r.post_refine->rotation_scale; - }; - constexpr int MAX_ROTATION_SCALE_ROUNDS = 8; - const RotationScaleWalk walk = WalkRotationScale(pass1.post_refine->rotation_scale, index_at, - refit_at, MAX_ROTATION_SCALE_ROUNDS); - experiment_ = before_walk; - prepass_detector_geometry_ = geometry_before_walk; - prepass_mosaicity_ = mosaicity_before_walk; - prepass_rotation_scale_.reset(); - if (!walk.trail.empty()) { - rotation_scale_walk = fmt::format( - "goniometer rotation scale walked from the pass-1 fit {:.5f}, validation spots on " - "the lattice beyond chance by scale: {} - {}", pass1.post_refine->rotation_scale, - walk.trail, walk.scale == 1.0f ? "the stored angles stand" - : fmt::format("{:.5f} adopted", walk.scale)); - logger.Info("Two-pass: {}", rotation_scale_walk); - } - note(walk.scale == 1.0f ? std::string("the stored rotation angles stand") - : fmt::format("rotation scale {:.5f} adopted", walk.scale)); - if (walk.scale != 1.0f) { - prepass_rotation_scale_ = walk.scale; - pass1.post_refine->rotation_scale = walk.scale; - pass1.post_refine->rotation_scale_suspect = true; - logger.Warning("Goniometer rotation scale {:.5f} ({:+.2f} %): the stage turned {} than the " - "angles stored in the file, which are the COMMANDED values. The second pass " - "integrates at the corrected angles; the fault is in the hardware and should " - "be fixed there", walk.scale, 100.0 * (walk.scale - 1.0), - walk.scale > 1.0f ? "further" : "less far"); - } - } - - // ... and apply the adopted scale, on the same measure-then-re-integrate footing as the geometry. - if (prepass_rotation_scale_ && gonio_snapshot) { - experiment_.Goniometer(ScaleRotation(*experiment_.GetGoniometer(), *prepass_rotation_scale_)); - // The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has - // just stopped using, and the override can only ever raise the second pass's own estimate (it - // takes the larger of the two). Carrying it over would hold the second pass at the rocking - // width the uncorrected angles produced - the correction half-applied. Drop it and let the - // second pass fit its own. - prepass_mosaicity_.clear(); - logger.Info("Two-pass: goniometer rotation scale {:.5f} will drive the second integration pass " - "(oscillation {:.4f} -> {:.4f} deg per image)", *prepass_rotation_scale_, - gonio_snapshot->GetIncrement_deg(), - gonio_snapshot->GetIncrement_deg() * *prepass_rotation_scale_); - } - // The integration radius the pre-scan measured from the crystal's own spot width, applied now - // that the geometry is fitted. The pre-pass ran at the radius the run started with (see - // bragg_adaptive_), so everything pass 2 consumes from it - the post-refined distance and beam, - // the space group, the smoothed mosaicity - is what the fixed radius produces, and the widened - // disk changes only the intensities of the canonical pass. - if (bragg_adaptive_) { - experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_); - logger.Info("Two-pass: the canonical pass integrates at the measured spot width, " - "r1={:.1f} r2={:.1f} r3={:.2f}", bragg_adaptive_->GetR1(), - bragg_adaptive_->GetR2(), bragg_adaptive_->GetR3()); - } - - // Let every refined pass from here on ask the quality guard's question of itself, so one that - // is going to be thrown away can skip the report and the files it would have superseded. - quality_guard_pass1_ = pass1; - config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix) - // The canonical pass measures the post-refinement again where the pre-pass moved the geometry - // (it applies nothing), so that the walk below has something to decide on. - postrefine_probe_ = prepass_detector_geometry_.has_value(); - - const auto current_geometry = [&] { - return std::array{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), - experiment_.GetDetectorDistance_mm(), - experiment_.GetPoniRot1_rad(), experiment_.GetPoniRot2_rad()}; - }; - const auto set_geometry = [&](const std::array &g) { - experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) - .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); - }; - - // The metric symmetry, asked of the spots rather than of the tolerance that admitted it. The - // Bravais class is chosen by a walk that reads two axes as EQUAL when they agree to a fixed - // relative tolerance, and the class carrying that equality is then imposed everywhere below: - // the cell is refined with a = b, reflections are predicted from it and integrated at the - // predicted positions. On a small cell the tolerance is far wider than the positions resolve - - // a 5.8 A axis pair 2.4 % apart is read as tetragonal, and the mean the constraint substitutes - // mispredicts the outermost reflections by a good ten pixels against an integration disc of - // four - so a genuinely orthorhombic crystal is integrated as tetragonal and every decision - // downstream, up to the space group, is read off the wreckage. - // - // So the equality is measured rather than assumed. The indexer already refines the same - // candidate a second time with nothing held, and pass 1 reports how far that free refinement - // leaves the two axes apart (lattice_length_split). A relative split of eps displaces a - // reflection at radius r by eps/2 * r pixels, so the question "do these data resolve it" has - // an answer in the units the integration works in: where that displacement at the far corner - // of the detector stays inside the integration disc, imposing the equality moves nothing out - // of its own box and the higher symmetry is kept with no pass run. Where it does not, both - // hypotheses are run - the promoted class, and the class the same walk carries when it is - // granted no length equality at all - and whichever REALISES the lower held-out residual on a - // re-integrated, re-indexed pass is kept, the - // promoted class on a tie, because it is the one that asserts more. Both arms are probe-only - // passes (postrefine_probe_only_): the residual is all the comparison reads and it is in hand - // before the scaling engine is built, so neither arm pays for a merge, a space-group search or - // a set of correction surfaces, and only the class that wins is ever run as a canonical pass. - std::string lattice_arm; // the arm decision, for the pass decision text - // - // An ANGLE promotion is the same question: the walk holds a monoclinic or orthorhombic angle at - // 90 when the unrefined cell puts it within LATTICE_SEARCH_ANGLE_TOLERANCE_DEG of it, and an - // angle error of delta displaces a reflection at radius r by up to delta/2 * r pixels. Its - // demoted arm is the indexer's free refinement itself; the space group, and with it any - // lower class the data do carry, is then found from the intensities as for any triclinic - // indexing. Measured: a monoclinic crystal whose beta sits 1.4 deg from 90 indexed every - // validation frame as primitive orthorhombic, so the per-frame guard never fired. - const bool length_arm = pass1.lattice_length_split.has_value(); - const bool angle_arm = !length_arm && pass1.lattice_angle_departure_deg.has_value(); - if (!cancelled_ && (length_arm || angle_arm) && pass1.rotation_lattice_type) { - const double split = length_arm ? *pass1.lattice_length_split - : *pass1.lattice_angle_departure_deg * PI / 180.0; - const double dx = std::max(experiment_.GetBeamX_pxl(), - static_cast(experiment_.GetXPixelsNum()) - experiment_.GetBeamX_pxl()); - const double dy = std::max(experiment_.GetBeamY_pxl(), - static_cast(experiment_.GetYPixelsNum()) - experiment_.GetBeamY_pxl()); - const double r_corner_px = std::hypot(dx, dy); - const double miss_px = 0.5 * split * r_corner_px; - const double disc_px = experiment_.GetBraggIntegrationSettings().GetR1(); - if (miss_px > disc_px) { - if (length_arm) - logger.Info("Metric symmetry: {} was promoted on a = b, but the free refinement leaves " - "the two axes {:.2f} % apart - {:.1f} px of misprediction at the far corner " - "against an integration disc of {:.1f} px. Both classes are re-integrated " - "and judged on the held-out residual each realises", - gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), - 100.0 * split, miss_px, disc_px); - else - logger.Info("Metric symmetry: {} holds its angles at 90 deg, but the free refinement " - "leaves one {:.2f} deg from it - {:.1f} px of misprediction at the far " - "corner against an integration disc of {:.1f} px. Both classes are " - "re-integrated and judged on the held-out residual each realises", - gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), - *pass1.lattice_angle_departure_deg, miss_px, disc_px); - step("Metric-symmetry check"); - note(fmt::format("{} re-integrated with and without its {} equality", - gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), - length_arm ? "a = b" : "90 deg angle")); - // The arms are compared on the post-refinement's held-out residual, so they have to - // measure one: a canonical pass only does where the pre-pass moved the geometry. - const bool probe_was = postrefine_probe_; - postrefine_probe_ = true; - postrefine_probe_only_ = true; - const auto promoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); - ++arm_passes; - IndexingSettings demoted_settings = experiment_.GetIndexingSettings(); - demoted_settings.LatticeClassDemoted(length_arm); - experiment_.ImportIndexingSettings(demoted_settings); - lattice_angles_demoted_ = angle_arm; - const auto demoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); - ++arm_passes; - postrefine_probe_only_ = false; - postrefine_probe_ = probe_was; - const double realised_promoted = promoted_arm.post_refine ? promoted_arm.post_refine->held_out_before : NAN; - const double realised_demoted = demoted_arm.post_refine ? demoted_arm.post_refine->held_out_before : NAN; - const auto class_text = [](const ProcessResult &r) { - return r.rotation_lattice_type - ? fmt::format("{}-centred {}", r.rotation_lattice_type->centering, - gemmi::crystal_system_str(r.rotation_lattice_type->crystal_system)) - : std::string("no lattice"); - }; - if (!cancelled_ && promoted_arm.post_refine && demoted_arm.post_refine - && HeldOutResidualFell(*promoted_arm.post_refine, *demoted_arm.post_refine)) { - logger.Warning("Metric symmetry: re-integrated as {} the held-out residual is {:.3e} " - "against {:.3e} as {} - the equality is not this crystal's, and the " - "run goes on without it", class_text(demoted_arm), realised_demoted, - realised_promoted, class_text(promoted_arm)); - note(fmt::format("the equality dropped: {} fits better", class_text(demoted_arm))); - lattice_arm = fmt::format( - "the metric symmetry was promoted to {} on {} the free refinement leaves " - "{:.2f} {}; re-integrated as {} the held-out residual is {:.3e} against " - "{:.3e}, so the equality was dropped", class_text(promoted_arm), - length_arm ? "axes" : "an angle", length_arm ? 100.0 * split : split * 180.0 / PI, - length_arm ? "% apart" : "deg from 90", - class_text(demoted_arm), realised_demoted, realised_promoted); - } else if (!cancelled_) { - note(fmt::format("the metric symmetry stands: {}", class_text(promoted_arm))); - IndexingSettings promoted_settings = experiment_.GetIndexingSettings(); - promoted_settings.LatticeClassDemoted(false); - experiment_.ImportIndexingSettings(promoted_settings); - lattice_angles_demoted_ = false; - // An arm that indexed nothing, or post-refined nothing, measures no residual. The - // metric symmetry stands either way - that is the arm whose absence of evidence is - // not evidence against it - but the run must not report a comparison it could not - // make: written as a number, a missing residual came out as "realised a held-out - // residual of nan against 1.234e-03, so the metric symmetry stands". - const bool compared = std::isfinite(realised_promoted) && std::isfinite(realised_demoted); - if (compared) { - logger.Info("Metric symmetry: re-integrated as {} the held-out residual is {:.3e}, " - "which {} does not beat ({:.3e}) - the higher symmetry stands", - class_text(promoted_arm), realised_promoted, class_text(demoted_arm), - realised_demoted); - lattice_arm = fmt::format( - "{} was run as well and realised a held-out residual of {:.3e} against {:.3e} " - "as {}, so the metric symmetry stands", class_text(demoted_arm), - realised_demoted, realised_promoted, class_text(promoted_arm)); - } else { - const char *which = !std::isfinite(realised_demoted) && !std::isfinite(realised_promoted) - ? "neither arm measured a held-out residual" - : !std::isfinite(realised_demoted) - ? "the demoted arm measured no held-out residual" - : "the promoted arm measured no held-out residual"; - logger.Info("Metric symmetry: {}, so the two could not be compared - the higher " - "symmetry stands", which); - lattice_arm = fmt::format( - "{} was run as well, but {}, so the two could not be compared and the " - "metric symmetry stands untested", class_text(demoted_arm), which); - } - } - } else { - logger.Info("Metric symmetry: {} - the free refinement departs from it by {:.2f} {}, " - "{:.1f} px at the far corner, inside the integration disc of {:.1f} px, so " - "these data do not resolve the difference", - gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), - length_arm ? 100.0 * split : split * 180.0 / PI, length_arm ? "%" : "deg", - miss_px, disc_px); - } - } - - logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*", - base_prefix); - speculate_geometry_probe_ = true; // a walk that follows starts with an indexing probe - rerun_if_starved_ = true; - auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); - rerun_if_starved_ = false; - speculate_geometry_probe_ = false; - - // Walk the geometry to the fit's fixed point. The post-refinement is fitted to reflections - // integrated at the very geometry it is correcting, so one fit reaches only as far as those - // reflections see: on a header whose distance is several percent long, the pre-pass absorbs the - // error into an inflated cell, and from that self-consistent pair the fit finds a fraction of it - // and stops. The canonical pass has just re-integrated and re-indexed at what it committed, so - // asking the same fit there is the next step of the same walk - and that re-indexing is also the - // RATIFICATION the old size bound stood in for: a second lattice does not index better at a - // geometry it did not come from, a wrong header does, and a pass that lost the crystal is thrown - // out by the quality guard below, which returns the run to the header geometry. - // - // What starts the walk outright is a move of more than one step: that is the fit saying the - // geometry is somewhere else entirely rather than a fraction of a percent away. A smaller move - // has to be preferred by the validation frames first (below), so a run whose fit settles beside - // the geometry in hand costs two indexing probes and no pass. Once the run IS walking - // it keeps walking because the steps do not get smaller in proportion to what is left: measured - // on a header 3.9 % long, the first round took a third of the error and each of the next took a - // third of the rest. - // - // But a round is kept only for what it REALISES, not for what the fit predicts: the held-out - // residual of the pass re-integrated and re-indexed at the new geometry has to fall below the - // best round's by more than the cross-validation noise of the two (HeldOutResidualFell). With - // the walk allowed twelve rounds, every round's fit predicted a gain the re-integration did not - // show - the realised residual stayed flat, or passed a minimum and rose again, while the - // distance crept on along the distance/cell-scale degeneracy - so "the fit still commits" never - // ended a walk, and the cap decided where on the drift the run stopped. A round that does not - // pay ends the walk and the run goes back to the best round it reached. So does a fit whose - // committed cell the re-indexing pushes back the same way twice running (ReindexPushesCellBack): - // the fit keeps asking for a cell the indexer at the fit's own geometry does not return, and - // nothing ratifies it. A walk still paying when the rounds run out keeps its best round - the - // last - and says it did not converge; the quality guard below judges it against the header - // pass as it judges every refined pass. - // - // The residual is not the only thing a round realises. It is the mean over every reflection, - // measured by centroids taken inside a disc centred on the prediction - so they follow the - // prediction part of the way - and dominated by the low-resolution reflections, where a distance - // and the cell scale that compensates it move every spot alike. What does see such an error is - // the re-indexing: the share of the validation spots that land on the lattice, which the - // high-resolution spots lose first. Measured on a crystal whose canonical pass ran 0.6 % long, - // the round at the fit's distance put 70.4 % of the validation spots on the lattice against - // 58.1 %, while reading the residual only 0.74 of its noise lower - and the run the residual - // stopped kept the long distance and lost 0.05 A of resolution. So a round is kept when EITHER - // realises a gain beyond its noise: the residual (HeldOutResidualFell), or the validation - // frames (ValidationEvidencePrefers). - // - // And a move of less than one step no longer settles the geometry by being small: it is asked - // of the validation frames too, but cheaply first - an indexing probe at the fit's geometry - // against one at the geometry in hand, each scoring the lattice and stopping there - and only a - // move the probes prefer pays for a re-integrated round. A move of a step or more is run as a - // round directly, as it always was, and so is every round of the walk it starts. - const auto on_lattice_pct = [](const ValidationSpotEvidence &e) { - return 100.0 * static_cast(e.on_lattice - e.by_chance) - / static_cast(std::max(1, e.spots)); - }; - // An indexing probe at a geometry: the validation evidence of the lattice indexed there, and - // nothing else - the pass stops once the lattice is scored, and the experiment is put back. - const auto index_at_geometry = [&](const std::array &g) { - const DiffractionExperiment before_probe = experiment_; - const bool probe_was = postrefine_probe_; - set_geometry(g); - CollectSpeculativeProbe(); - indexing_probe_only_ = true; - postrefine_probe_ = false; - ProcessResult r; - try { - r = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Info("Two-pass: the indexing probe at distance {:.3f} mm did not complete ({})", - g[2], e.what()); - } - indexing_probe_only_ = false; - postrefine_probe_ = probe_was; - experiment_ = before_probe; - ++arm_passes; - return r.validation_evidence; - }; - constexpr int MAX_GEOMETRY_ROUNDS = 8; - int geometry_rounds = 0; // passes run at a geometry the walk moved to - int walk_passes = 0; // every pass the walk ran, the return to its best round included - int best_round = 0; - std::array best_geometry = current_geometry(); - PostRefineResult best_fit; // the fit measured in the best round's pass: its realised residual - int last_pushback = 0; - bool geometry_not_converged = false; - std::string walk_stop; // why the walk ended, where anything but the fit settling ended it - std::string walk_rounds; // distance / realised residual per round, for the pass decision - if (pass2.post_refine) { - best_fit = *pass2.post_refine; - walk_rounds = fmt::format("{:.3f} mm / {:.3e} / {:.1f}%", best_geometry[2], best_fit.held_out_before, - on_lattice_pct(pass2.validation_evidence)); - } - ValidationSpotEvidence best_evidence = pass2.validation_evidence; - std::optional probe_at_best; // the indexing probe at best_geometry, once run - bool large_walk = false; // a round of this walk was started by a move of a step or more - // What the last round started from: a pass can leave the run changed behind it (a spot budget - // it shortened, an axis sign it flipped), so a round is written - // from here, not from where its own probe left the run. - DiffractionExperiment last_round_start = experiment_; - bool walk_shown = false; - while (!cancelled_ && prepass_detector_geometry_ - && pass2.post_refine && pass2.post_refine->detector_refined) { - if (!walk_shown) { - step("Geometry check"); - walk_shown = true; - } - if (geometry_rounds == MAX_GEOMETRY_ROUNDS) { - geometry_not_converged = true; - walk_stop = fmt::format("the fit still moved the geometry after {} rounds, every one of " - "which lowered the realised residual", MAX_GEOMETRY_ROUNDS); - logger.Warning("Two-pass: {} - the walk did not converge; the run keeps its last round " - "and the quality guard judges it against the header pass", walk_stop); - break; - } - const PostRefineResult fit = *pass2.post_refine; - const std::array g = *prepass_detector_geometry_; - large_walk = large_walk || fit.large_move; - if (!large_walk) { - if (!probe_at_best) - probe_at_best = index_at_geometry(best_geometry); - const ValidationSpotEvidence probe_at_fit = index_at_geometry(g); - const bool prefers = ValidationEvidencePrefers(*probe_at_best, probe_at_fit); - logger.Info("Two-pass: the post-refinement moves the geometry by less than a step (distance " - "{:.3f} -> {:.3f} mm, beam {:.2f},{:.2f} -> {:.2f},{:.2f} px); indexing probes put " - "{:.1f}% of the validation spots on the lattice beyond chance there against {:.1f}% " - "at the geometry in hand - {}", fit.distance_before_mm, fit.distance_after_mm, - fit.beam_x_before_px, fit.beam_y_before_px, fit.beam_x_after_px, fit.beam_y_after_px, - on_lattice_pct(probe_at_fit), on_lattice_pct(*probe_at_best), - prefers ? "worth a round" : "the geometry in hand stands"); - if (!prefers) - break; - probe_at_best = probe_at_fit; - } - logger.Info("Two-pass: the post-refinement at the adopted geometry still moves it " - "(distance {:.3f} -> {:.3f} mm) - re-integrating and re-indexing at what it " - "asks for (round {} of at most {})", fit.distance_before_mm, - fit.distance_after_mm, geometry_rounds + 1, MAX_GEOMETRY_ROUNDS); - set_geometry(g); - last_round_start = experiment_; - // A round is a probe pass: the walk reads only its fit and its validation evidence, both in - // hand before the merge, so the files are written once, below, at the round the walk keeps. - speculate_geometry_probe_ = !large_walk; - postrefine_probe_only_ = true; - pass2 = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); - postrefine_probe_only_ = false; - speculate_geometry_probe_ = false; - ++geometry_rounds; - ++walk_passes; - const double realised = pass2.post_refine ? pass2.post_refine->held_out_before : NAN; - walk_rounds += fmt::format(", {:.3f} mm / {:.3e} / {:.1f}%", g[2], realised, - on_lattice_pct(pass2.validation_evidence)); - const bool residual_fell = pass2.post_refine && HeldOutResidualFell(best_fit, *pass2.post_refine); - const bool evidence_prefers = ValidationEvidencePrefers(best_evidence, pass2.validation_evidence); - if (!pass2.post_refine || !(residual_fell || evidence_prefers)) { - walk_stop = fmt::format( - "round {} lowered neither the realised held-out residual by more than its noise " - "({:.3e} against {:.3e}, standard errors {:.1e} and {:.1e}) nor the share of " - "validation spots off the lattice ({:.1f}% on it beyond chance against {:.1f}%)", - geometry_rounds, realised, best_fit.held_out_before, - pass2.post_refine ? pass2.post_refine->held_out_before_se : NAN, - best_fit.held_out_before_se, on_lattice_pct(pass2.validation_evidence), - on_lattice_pct(best_evidence)); - logger.Info("Two-pass: {}", walk_stop); - break; - } - best_round = geometry_rounds; - best_geometry = g; - best_fit = *pass2.post_refine; - best_evidence = pass2.validation_evidence; - const int pushback = ReindexPushesCellBack(fit, *pass2.post_refine); - if (pushback != 0 && pushback == last_pushback) { - walk_stop = fmt::format( - "re-indexing pushed the cell length {} back against the fit for the second " - "round running, so nothing ratifies the cell the fit asks for", - "abc"[std::abs(pushback) - 1]); - logger.Info("Two-pass: {} - the walk stops at round {}", walk_stop, geometry_rounds); - break; - } - last_pushback = pushback; - } - postrefine_probe_ = false; - if (walk_shown) - note(geometry_rounds == 0 ? std::string("the geometry in hand stands") - : fmt::format("geometry walked {} round{}, round {} kept: distance {:.2f} mm{}", geometry_rounds, - geometry_rounds == 1 ? "" : "s", best_round, best_geometry[2], - best_round > 0 ? " - the files are written again there" : "")); - // The rounds wrote nothing, so the round the walk kept is run again to write the canonical files. - if (!cancelled_ && geometry_rounds > 0) { - if (best_round == geometry_rounds) { - logger.Info("Two-pass: writing round {} of the geometry walk (distance {:.3f} mm)", - best_round, best_geometry[2]); - experiment_ = last_round_start; - } else { - logger.Info("Two-pass: going back to round {} of the geometry walk (distance {:.3f} mm), " - "the last one that lowered the realised residual", best_round, best_geometry[2]); - set_geometry(best_geometry); - } - rerun_if_starved_ = true; - pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); - rerun_if_starved_ = false; - ++walk_passes; - } - - // The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever - // result is returned - it is a result of the run, not of the pass that happened to fit it. Where - // the walk took rounds they are ONE geometry decision, so what the report names is where the - // geometry started and where the run left it, not the first step of the way there. - if (best_round > 0 && pass1.post_refine) { - pass1.post_refine->beam_x_after_px = best_geometry[0]; - pass1.post_refine->beam_y_after_px = best_geometry[1]; - pass1.post_refine->distance_after_mm = best_geometry[2]; - // A first pass that held the header distance and a walk that then left it: the run's - // decision is the walk's. - pass1.post_refine->distance_held = false; - } - pass2.post_refine = pass1.post_refine; - pass2.pass_number = 2 + arm_passes + walk_passes; - pass2.pass_count = 2 + arm_passes + walk_passes; - if (geometry_rounds > 0) - pass2.pass_decision = fmt::format( - "post-refined geometry walked {} round{} and kept round {}, detector distance {:.3f} mm " - "against the header's {:.3f}: {}; distance / realised held-out residual / validation spots on the lattice beyond chance by round: {}", - geometry_rounds, geometry_rounds == 1 ? "" : "s", best_round, best_geometry[2], - pass1.post_refine ? pass1.post_refine->distance_before_mm : 0.0, - walk_stop.empty() ? "the fit at the kept geometry commits no further move" : walk_stop, - walk_rounds); - - // r2 = r1 + 2 is also the inner edge of the background ring, so widening the signal disk pushes - // that ring out into the neighbours; a reflection whose ring is left with five or fewer clean - // pixels has no background and is dropped whole. On a pattern dense enough in three dimensions - // at once that costs a quarter of the observations - and NOTHING available before integrating - // predicts it: the crystal it happens to has a mid-table predicted spot spacing, wider than - // crystals that survive r1 = 12. What does measure it is the integrator counting its own drops. - // The pass that integrates at the widened radius is now the canonical one, so the measurement - // and the remedy are both here: re-run it at the radius that keeps the data, the same remedy the - // two guards below use and, like them, a pass only on a crystal that was going to be wrong. - // Asked of every pass that can end the run - the two guards below re-run at whatever radius is - // in hand, and a starved refined pass has written nothing (see the superseded flag in - // RunPipeline) - so each is followed by it. After the first re-run the radius is the fixed one - // and this never fires again. - const auto rerun_if_starved = [&](ProcessResult &pass, const std::string decision) { - if (cancelled_ || !bragg_adaptive_ || !bkg_starved_fraction_ - || *bkg_starved_fraction_ <= spot_width::BKG_STARVED_MAX_FRACTION) - return; - logger.Info("Two-pass: at the widened integration radius r1={:.1f} the neighbouring " - "reflections left {:.2f}% of the predicted reflections without a background " - "ring (bound {:.2f}%), so this pass is re-run at r1={:.1f} r2={:.1f} r3={:.2f}", - bragg_adaptive_->GetR1(), 100.0 * *bkg_starved_fraction_, - 100.0 * spot_width::BKG_STARVED_MAX_FRACTION, - bragg_before_adaptive_->GetR1(), bragg_before_adaptive_->GetR2(), - bragg_before_adaptive_->GetR3()); - const double starved = *bkg_starved_fraction_; - const float wide_r1 = bragg_adaptive_->GetR1(); - const int passes = pass.pass_count + 1; - experiment_.ImportBraggIntegrationSettings(*bragg_before_adaptive_); - bragg_adaptive_.reset(); - note(fmt::format("{:.1f}% of the background rings starved: the pass runs again at r1 {:.1f} px", - 100.0 * starved, bragg_before_adaptive_->GetR1())); - pass = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); - pass.post_refine = pass1.post_refine; - pass.pass_number = passes; - pass.pass_count = passes; - pass.pass_decision = fmt::format( - "{}, re-run at the fixed integration radius: at r1={:.1f} neighbouring reflections " - "left {:.2f}% of the predicted reflections without a background ring", decision, wide_r1, - 100.0 * starved); - }; - rerun_if_starved(pass2, "post-refined geometry adopted"); - - // Supercell-collapse guard for the second pass. The de-novo lattice search can be BISTABLE on a - // periodically-modulated dataset: the tiny post-refine geometry move tips it out of the true-cell - // basin onto a spurious SUPERCELL (a large multiple of pass-1's cell volume), and the space group - // reused from pass-1 is then stamped onto that wrong-shaped cell - catastrophic. Pass 1 (header - // geometry, un-perturbed) found the correct lattice, so re-run the second pass with pass-1's FULL - // indexing result forced (its refined lattice + orientation + rotation axis + geometry - not just the - // cell, which loses the per-frame orientation and integrates worse): the re-run then matches de-novo - // quality but at the refined geometry. Only fires when pass 2 is markedly LARGER than pass 1 - a - // genuine de-novo demotion to a smaller primitive (the pseudo-symmetry case the second pass re-indexes - // de novo to catch) has pass2 <= pass1 and is left untouched. - if (!cancelled_ && prepass_result_ && pass1.consensus_cell && pass2.consensus_cell) { - // PRIMITIVE volumes, like the two guards in RunPipeline: a centred conventional cell is an - // exact integer multiple of its primitive one (C/I 2x, R 3x, F 4x), so two settings of the - // SAME lattice differ by that factor and comparing the conventional cells reads a mere - // change of setting as a supercell - which then forces pass 1's setting and, with it, its - // lower symmetry. Measured: every dataset this guard fired on was at an exact centring - // multiplicity, one of them an F-centred lattice held down to P1 at 3.99x. - const auto primitive_volume = [](const UnitCell &cell, std::optional centering) { - return std::abs(CrystalLattice(cell).ToPrimitive(centering.value_or('P')).CalcVolume()); - }; - const double v1 = primitive_volume(*pass1.consensus_cell, pass1.consensus_centering); - const double v2 = primitive_volume(*pass2.consensus_cell, pass2.consensus_centering); - if (v1 > 0.0 && v2 > 1.5 * v1) { - logger.Info("Two-pass: second-pass primitive cell volume {:.0f} A^3 is {:.2f}x pass-1's " - "({:.0f} A^3) - a spurious supercell; re-running the second pass with " - "pass-1's result forced", - v2, v2 / v1, v1); - note("a spurious supercell: the pass runs again on the pre-pass lattice"); - force_rotation_result_ = *prepass_result_; - // Pass 1's result carries pass 1's goniometer, and the per-image path takes its angles - // from there - forcing it whole would put the uncorrected angles back and silently undo - // the rotation scale for the re-run. - if (prepass_rotation_scale_ && force_rotation_result_->axis) - force_rotation_result_->axis = - ScaleRotation(*force_rotation_result_->axis, *prepass_rotation_scale_); - rerun_if_starved_ = true; - pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); - rerun_if_starved_ = false; - pass2.post_refine = pass1.post_refine; - pass2.pass_number = 3; - pass2.pass_count = 3; - pass2.pass_decision = fmt::format( - "post-refined geometry adopted, re-run with pass-1's lattice forced: the de-novo " - "primitive cell volume {:.0f} A^3 was {:.2f}x pass-1's ({:.0f} A^3), a spurious " - "supercell", v2, v2 / v1, v1); - // Still forced: the re-run at the fixed radius is the same pass. - rerun_if_starved(pass2, pass2.pass_decision); - force_rotation_result_.reset(); - } - } - // The refined pass is normally the better answer, which is why it is the canonical output - - // but it is not guaranteed to be, and the guard that decides is RefinedPassIsWorse above. It - // has already been asked once, inside the pass itself, which is how a pass that is about to be - // thrown away knows not to write a report; this is where the answer is acted on. - if (!cancelled_ && pass1.has_merge_statistics && pass2.has_merge_statistics) { - logger.Info("Two-pass quality guard: the refined pass holds {} reflections at I/sigma >= 2 and " - "{} at <= -2 of {} in a {:.0f} A^3 primitive cell, the header-geometry pass {} and {} " - "of {} in {:.0f} A^3; low-order axial {} against {}", - pass2.search_merge_strong_reflections, pass2.search_merge_negative_reflections, - pass2.search_merge_reflections, PrimitiveCellVolume(pass2), - pass1.search_merge_strong_reflections, pass1.search_merge_negative_reflections, - pass1.search_merge_reflections, PrimitiveCellVolume(pass1), - pass2.search_merge_axial_reflections, pass1.search_merge_axial_reflections); - const std::string worse = RefinedPassIsWorse(pass1, pass2); - if (!worse.empty()) { - logger.Warning("Two-pass: the refined pass is worse than the header-geometry pass ({}) " - "- going back to the header geometry. The refined geometry did not help " - "this crystal.", worse); - // The refined pass has already written the canonical files, so re-run at the header - // geometry to replace them - the same remedy the supercell collapse above uses, and it - // only costs a pass on a crystal that was going to be wrong otherwise. - experiment_.BeamX_pxl(header_geometry[0]).BeamY_pxl(header_geometry[1]) - .DetectorDistance_mm(header_geometry[2]) - .PoniRot1_rad(header_geometry[3]).PoniRot2_rad(header_geometry[4]); - if (prepass_rotation_scale_ && gonio_snapshot) { - // The header rotation angles with it - but not the axis DIRECTION, which is not - // part of the geometry being reverted and may have been corrected by the sign - // rescue (see the restore before pass 2). - GoniometerAxis restored = *gonio_snapshot; - if (const auto current = experiment_.GetGoniometer()) - restored.Axis(current->GetAxis()); - experiment_.Goniometer(restored); - prepass_rotation_scale_.reset(); - pass1.post_refine->rotation_scale_suspect = false; - } - config_.output_prefix = base_prefix; - // This pass is the answer whatever it measures - the guard has had its one chance - - // so it must not conclude from its own search merge that it is going to be re-run. - quality_guard_pass1_.reset(); - // What the guard preferred is pass 1 - its lattice, orientation and axis at this very - // geometry - so that is what is re-integrated, not a fresh de-novo index. A de-novo - // re-index here is a new hypothesis nobody judged: measured, at the same centre where - // pass 1 indexed 56/60 validation frames it found 5/60, took the opposite axis sign and - // shipped a merge with R_meas in the hundreds of percent. - if (prepass_result_) - force_rotation_result_ = *prepass_result_; - rerun_if_starved_ = true; - note("the refined geometry is worse: the pass runs again at the header geometry"); - auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); - rerun_if_starved_ = false; - redo.post_refine = pass1.post_refine; - redo.pass_number = pass2.pass_count + 1; - redo.pass_count = pass2.pass_count + 1; - redo.pass_decision = fmt::format( - "header geometry re-adopted: the post-refined pass was worse ({})", worse); - // Still forced: the re-run at the fixed radius is the same pass. - rerun_if_starved(redo, redo.pass_decision); - force_rotation_result_.reset(); - return redo; - } - if (pass2.pass_decision.empty()) - pass2.pass_decision = fmt::format( - "post-refined geometry adopted ({}, {})", SearchMergeCompletenessText(pass1, pass2), - SignalText(pass1, pass2)); - } - if (pass2.pass_decision.empty()) - pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_) - ? "post-refined geometry adopted" - : "the post-refinement committed no geometry change, so this pass reproduces the first"; - if (!lattice_arm.empty()) - pass2.pass_decision = lattice_arm + "; " + pass2.pass_decision; - if (!rotation_scale_walk.empty()) - pass2.pass_decision = rotation_scale_walk + "; " + pass2.pass_decision; - pass2.geometry_not_converged = geometry_not_converged; - return pass2; - } - return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); -} - namespace { // Angle between the crystal's symmetry axes and the spindle. A rotation sweep never records the // reflections whose reciprocal vector lies close to the spindle - the blind cusp - and normally diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp new file mode 100644 index 000000000..2bcd46a4c --- /dev/null +++ b/rugnux/RugnuxPasses.cpp @@ -0,0 +1,1242 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "Rugnux.h" +#include "../image_analysis/structure_refinement/ModelValidation.h" +#include "../writer/WriteModel.h" +#include "../image_analysis/indexing/SpindleCuspLoss.h" +#include "../image_analysis/bragg_integration/SpotWidth.h" +#include "../image_analysis/bragg_integration/SpotFootprint.h" +#include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA +#include "../image_analysis/hot_pixels/HotPixels.h" +#include "../image_analysis/SensorAbsorption.h" +#include "DiagnosticOutput.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../reader/JFJochHDF5Reader.h" +#include "../common/JFJochMath.h" +#include "../common/ParallelFor.h" +#include "../common/Logger.h" +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/AzimuthalIntegrationProfile.h" +#include "../common/CUDAWrapper.h" +#include "../common/JFJochException.h" +#include "../common/time_utc.h" +#include "../writer/FileWriter.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/beam_stop/ShadowFinder.h" +#include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" +#include "../image_analysis/indexing/HarmonicContamination.h" +#include "../image_analysis/geom_refinement/BeamCenterFFT.h" +#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" +#include "../image_analysis/geom_refinement/GeometryRefiner.h" +#include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" +#include "../image_analysis/spot_finding/SpotUtils.h" +#include "gemmi/twin.hpp" +#include "../image_analysis/azint/AzIntEngineCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderGPU.h" +#endif +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/RfreeFlags.h" +#include "../image_analysis/scale_merge/RotationScaleMerge.h" +#include "../image_analysis/scale_merge/ResolutionCutoff.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../image_analysis/scale_merge/ScalingResult.h" +#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/geom_refinement/PostRefine.h" +#include "../image_analysis/lattice_search/LatticeSearch.h" +#include "../image_analysis/lattice_search/LePageLattice.h" +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" +#include "../image_analysis/scale_merge/FrenchWilson.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" +#include "../image_analysis/scale_merge/TranslationalNCS.h" +#include "../image_analysis/scale_merge/HKLKey.h" +#include "../image_analysis/scale_merge/ScaleOnTheFly.h" +#include "../image_analysis/scale_merge/StillsPartialityRefine.h" +#include "../image_analysis/WriteReflections.h" +#include "../image_analysis/bragg_integration/CalcISigma.h" +#include "../common/Definitions.h" +#include "../common/CorrelationCoefficient.h" +#include +#include +#include "RugnuxPipeline.h" + +using namespace rugnux_internal; + +namespace { + // How much better the alternative beam centre's search merge has to be for the run to move off + // the centre the file gave it (the beam-centre arbitration in RunAllPasses). + constexpr double MAX_CC_HALF_LOSS = 0.05; + + // Say so rather than printing a zero that reads as a measurement. ("not compared", not "not + // measured": completeness IS measured and reported; what did not happen is its use in arbitrating + // between the two passes.) + std::string SearchMergeCompletenessText(const ProcessResult &pass1, const ProcessResult &pass2) { + if (!pass1.search_merge_completeness_measured || !pass2.search_merge_completeness_measured) + return "completeness not compared"; + return fmt::format("completeness {:.1f}% vs {:.1f}%", pass2.search_merge_completeness, + pass1.search_merge_completeness); + } + + // The volume of the pass's primitive cell, 0 where it has no cell. + double PrimitiveCellVolume(const ProcessResult &r) { + if (!r.consensus_cell) + return 0.0; + return std::abs(CrystalLattice(*r.consensus_cell).ToPrimitive(r.consensus_centering.value_or('P')) + .CalcVolume()); + } + + // The signal the two passes measured, as the quality guard reads it (RefinedPassIsWorse). + std::string SignalText(const ProcessResult &pass1, const ProcessResult &pass2) { + return fmt::format("{} reflections at I/sigma >= 2 less the {} at <= -2 before corrections, " + "against {} less {}", pass2.search_merge_strong_reflections, + pass2.search_merge_negative_reflections, pass1.search_merge_strong_reflections, + pass1.search_merge_negative_reflections); + } + + // Whether the two passes describe their lattice in the same setting, which is what makes a count + // of reflections at a given Miller index the same count on both sides. A pass that re-indexes at + // the refined geometry can settle on a centred setting of the same lattice the other pass held + // primitive; the cells are then different bases and only the counts taken over all reflections + // survive the change. + bool SameSetting(const ProcessResult &pass1, const ProcessResult &pass2) { + if (pass1.consensus_centering != pass2.consensus_centering + || !pass1.consensus_cell || !pass2.consensus_cell) + return false; + const UnitCell &c1 = *pass1.consensus_cell, &c2 = *pass2.consensus_cell; + const auto same_length = [](double a, double b) { return std::abs(a - b) <= 0.05 * (a + b) / 2; }; + return same_length(c1.a, c2.a) && same_length(c1.b, c2.b) && same_length(c1.c, c2.c); + } +} + +// Pass 2 is normally the better answer, which is why it is the canonical output - but it is not +// guaranteed to be. Returns why the refined pass is worse, in the words the run reports it in, or +// an empty string when it is not. +// +// The comparison is made on each pass's SEARCH merge, not on the final one. Each pass determines +// its own space group, so the two final merges can be in different groups, and a completeness or a +// correlation in different groups are not the same measurement. The search merge is P1 on both +// sides (or the group the user fixed), over the full resolution range and before any correction +// surface is fitted - and it is counted over the centring-allowed reflections on both sides, so +// that a pass which finds a centring the other did not is not judged on the extinguished half it +// merged as noise (see the capture in RunPipeline). +// +// What the two passes are compared on there is the SIGNAL each measured - the reflections it +// merged at I/sigma >= 2 - and not, as this guard decided until it was measured, that merge's +// CC1/2. Two measurements settled it: +// - Against an EXTERNAL arbiter (the site's own XDS merge of the same images, scored over a +// resolution range common to the arms) the count of unique reflections at I/sigma >= 2 agreed +// with which geometry is the more accurate one on 23 of 27 arm-dataset pairs. R_meas, CC1/2 +// and ISa - which measure whether a merge agrees with ITSELF - agreed on 13 of 27, a coin flip. +// - This merge's CC1/2 is pooled over the whole range, uncut and uncorrected, so the shells with +// no signal in them dominate it and they are exactly the shells a geometry move disturbs: it +// reads 0.13 on a crystal whose data merge at 0.995, and 0.93 against 0.76 on two runs whose +// merged data agree to 0.7% of R_meas. Deciding on it sent two crystals back to a geometry +// that ships 22% fewer reflections than the one they refused. +// +// Three ways the refined pass can be worse, then: it merges more unique reflections than the cell +// it settled on can hold, it measured decisively less signal, or it lost the axial rows. +// +// The completeness arm only bites when the search merge is in a fixed group (-S). De novo that +// merge is in P1 and does not count its possible reflections at all, so there is no completeness +// and it does not arm. Counting them there was measured and not taken: it costs 5.2 ms, but a +// merge that fills its own asymmetric unit reads 27-43% of the P1 hemisphere over 28 rotation +// crystals, so the 100.5% bound is never approached and no decision changes - it would only put a +// figure into the report that reads as the dataset's completeness and is not one. +// +// Asked twice, from the same numbers: in RunAllPasses, where the answer is acted on, and inside a +// refined pass as soon as its own search merge is in - a pass that already knows it will be thrown +// away leaves its report and its files unwritten (see the superseded flag in RunPipeline). +std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &pass2) { + constexpr double MAX_CREDIBLE_COMPLETENESS = 100.5; // rounding headroom, nothing more + // How much of the signal the refined pass has to lose before the run moves back off the + // geometry it measured. This is not a quality bar: it is set above what the two passes differ + // by for reasons that are not the geometry, and they never integrate quite the same way - the + // canonical pass integrates at the measured spot width where the pre-pass integrated at the + // default, and it predicts with the pre-pass's smoothed mosaicity. Measured over ~100 rotation + // crystals, that alone moves the count by -2% to +35%, and on a crystal where the two passes + // ran at the SAME geometry it moved it +10% (and the CC1/2 this guard used to decide on by + // -0.25). Below this, the difference says nothing about the geometry, and a comparison that + // cannot tell keeps the geometry the run measured. + constexpr double DECISIVE_SIGNAL_LOSS = 0.10; + // The signal is the count of reflections at I/sigma >= 2 LESS the count at I/sigma <= -2. An + // empty reflection - one the lattice predicts and the crystal does not diffract into - is as + // likely to land in either tail, so the negative tail is the merge's own measurement of how + // much of its strong count is noise, and what is left is what the crystal put there. It is what + // makes the two counts comparable when the passes hold different lattices: a pass on an n-fold + // supercell merges n times the reflections, n - 1 of every n of them empty, and their noise alone + // reached the strong count - measured, a header pass on a 7x supercell out-counted the refined + // pass on the crystal's own lattice and the run went back to the supercell, while that merge + // held more reflections at <= -2 than at >= 2. On a real superstructure the extra classes hold + // signal, and it is counted. + const auto signal = [](const ProcessResult &r) { + return static_cast(r.search_merge_strong_reflections - r.search_merge_negative_reflections); + }; + // ... and the loss has to be a loss rather than a dilution. The wider disk pulls weaker + // reflections into the merge and the error model is then refitted over them, so a pass can hold + // MORE reflections and fewer of them at I/sigma >= 2 without having measured less of the + // crystal: measured, one crystal merges 2.1x the reflections in its refined pass, a third of + // the strong ones, and its refined pass is the better answer by every final statistic. Where + // the refined pass covers no more of reciprocal space and measured less signal, the loss is a + // loss. Between two lattices the raw counts are no comparison: a lattice of n times the primitive + // volume holds n times the reflections to the same resolution, so against a header pass on a + // supercell the refined pass always held fewer and the clause never protected it. The count of + // the smaller lattice is scaled by the index n, which is an integer - the same lattice, whose + // volume two refinements put a percent or two apart, is index 1 and compares as it always did. + const auto coverage_does_not_grow = [&] { + double n1 = static_cast(pass1.search_merge_reflections); + double n2 = static_cast(pass2.search_merge_reflections); + const double v1 = PrimitiveCellVolume(pass1), v2 = PrimitiveCellVolume(pass2); + if (v1 > 0.0 && v2 > 0.0) { + if (v1 >= v2) + n2 *= std::round(v1 / v2); + else + n1 *= std::round(v2 / v1); + } + return n2 <= n1; + }; + + // The other arm: the axial rows. A bulk statistic is an average over tens of thousands of + // reflections and the principal axial rows are a few dozen of them, so a pass can measure half + // the row and move it by nothing - yet those are the reflections the systematic absences are + // read off, and losing them costs a screw axis and with it the space group. It happens because + // a reflection on an axial row crosses the Ewald sphere over a narrow range of the sweep, where + // a general reflection is spread over it: a pass that loses a wedge of the sweep loses whole + // axial rows while its bulk statistics improve. Measured over a dozen rotation crystals the two + // passes agree on this count to within 9% and usually exactly; the crystal this arm was written + // for measured 44%. + // + // It is asked only where the two counts are counts of the SAME reflections: the a and b of a + // C-centred setting are the face diagonals of the primitive cell, so "h00" is a different + // direction in the two merges and the two rows are not comparable at all. And it needs pass 1 + // to have measured enough of the rows for the difference to mean anything - below eight there + // is no zone a screw could be claimed on either side. A pass that measured decisively more + // signal has said something about the bulk of the data that a few dozen reflections cannot + // answer, and is kept. + constexpr double MIN_AXIAL_ROWS_KEPT = 0.75; + constexpr int64_t MIN_AXIAL_ROWS_TO_COMPARE = 8; + const double signal1 = signal(pass1); + const double signal2 = signal(pass2); + const int64_t axial1 = pass1.search_merge_axial_reflections; + const int64_t axial2 = pass2.search_merge_axial_reflections; + const bool measured_less = signal1 > 0 && signal2 < (1.0 - DECISIVE_SIGNAL_LOSS) * signal1 + && coverage_does_not_grow(); + const bool measured_more = signal1 > 0 && signal2 > (1.0 + DECISIVE_SIGNAL_LOSS) * signal1; + const bool lost_axial_rows = SameSetting(pass1, pass2) + && axial1 >= MIN_AXIAL_ROWS_TO_COMPARE + && axial2 < MIN_AXIAL_ROWS_KEPT * static_cast(axial1) + && !measured_more; + const bool worse = pass2.lattice_conflicts_with_fixed_sg + || (pass1.search_merge_completeness_measured + && pass2.search_merge_completeness_measured + && pass2.search_merge_completeness > MAX_CREDIBLE_COMPLETENESS) + || measured_less || lost_axial_rows; + if (!worse) + return {}; + return lost_axial_rows && !measured_less + ? fmt::format("it holds {} of the {} low-order axial reflections the header-geometry " + "pass measured, and its {} do not say it is better", axial2, axial1, + SignalText(pass1, pass2)) + : fmt::format("{}, {}", SearchMergeCompletenessText(pass1, pass2), SignalText(pass1, pass2)); +} + +bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b) { + if (!a.rotation_lattice_type || !b.rotation_lattice_type || !a.consensus_cell || !b.consensus_cell) + return false; + const auto primitive_volume = [](const ProcessResult &r) { + return std::abs(CrystalLattice(*r.consensus_cell) + .ToPrimitive(r.consensus_centering.value_or('P')).CalcVolume()); + }; + const double va = primitive_volume(a), vb = primitive_volume(b); + return a.rotation_lattice_type->crystal_system == b.rotation_lattice_type->crystal_system + && a.rotation_lattice_type->centering == b.rotation_lattice_type->centering + && std::max(va, vb) < 1.02 * std::min(va, vb); +} + +std::string MeasuredCentreWins(const ProcessResult &file_arm, const ProcessResult &measured_arm) { + if (!file_arm.has_merge_statistics || !measured_arm.has_merge_statistics) + return {}; + // Same lattice: which geometry is the more accurate, read off the signal each arm measured - the + // statistic the two-pass quality guard found to agree with an external arbiter where CC1/2 did not + // (RefinedPassIsWorse). Without its margin: that margin covers what the header and refined passes + // differ by for reasons that are not the geometry - a different spot width and mosaicity - and the + // two arms here integrate with the same settings, so the geometry is the only difference, and the + // file's centre has no standing the measured one lacks - neither indexed the whole spot list. + if (ArmsHoldSameLattice(file_arm, measured_arm)) { + if (measured_arm.search_merge_strong_reflections > file_arm.search_merge_strong_reflections) + return fmt::format("the measured centre measures more of the crystal - {} reflections at " + "I/sigma >= 2 before corrections against {}", + measured_arm.search_merge_strong_reflections, + file_arm.search_merge_strong_reflections); + // ...or its strongest reflections agree with themselves where the file's do not. The pooled + // CC1/2 is no arbiter here - the shells with no signal in them dominate it - but in the + // lowest-resolution shell every reflection is strong, and half-sets that disagree there are + // equivalents given different indices: a centre about one reflection spacing out. Measured, + // such a file centre merged 2 % MORE reflections at I/sigma >= 2 than the right one, at an + // inner-shell CC1/2 of 0.610 against 0.997. + if (measured_arm.search_merge_inner_cc_half > file_arm.search_merge_inner_cc_half + MAX_CC_HALF_LOSS) + return fmt::format("the measured centre's strongest reflections agree where the file's do not - " + "CC1/2 {:.3f} against {:.3f} in the lowest-resolution shell of the search merge", + measured_arm.search_merge_inner_cc_half, file_arm.search_merge_inner_cc_half); + return {}; + } + if (measured_arm.search_merge_cc_half > file_arm.search_merge_cc_half + MAX_CC_HALF_LOSS) + return fmt::format("the measured centre merges better - CC1/2 before corrections {:.3f} against " + "{:.3f} on the search merge", measured_arm.search_merge_cc_half, + file_arm.search_merge_cc_half); + return {}; +} + +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. + quality_guard_pass1_.reset(); + // Rotation two-pass geometry post-refinement: the first pass integrates at the header geometry and + // post-refines the detector geometry (distance + beam from the observed spot positions, cell scale + axis + // from phi_obs); the second pass re-indexes and re-integrates with the refined geometry. Only the second + // (refined-geometry) pass writes merged files, under the plain "_*" name. The first pass still runs + // under a "_01" prefix so that the _process.h5 it writes, when one is asked for, does not collide + // with the second pass's - the ice-ring flags and geometry stored in each are the ones that pass used. + if (config_.mode == ProcessMode::FullAnalysis && config_.rotation_postrefine_geometry + && experiment_.IsRotationIndexing()) { + Logger logger("Rugnux"); + // The checks run between the passes are steps of their own on the console; their probe passes + // post none. + const auto step = [observer](const std::string &title) { if (observer) observer->OnStep(title); }; + const auto note = [observer](const std::string &text) { if (observer) observer->OnNote(text); }; + const std::string base_prefix = config_.output_prefix; + const auto gonio_snapshot = experiment_.GetGoniometer(); + // What the first pass starts from, so a second first pass at a different beam centre can start + // from the same place (the beam-centre arbitration below). + const DiffractionExperiment experiment_before_first_pass = experiment_; + const int64_t max_spot_count_snapshot = experiment_.GetDatasetSettings().GetMaxSpotCount(); + header_distance_mm_ = experiment_.GetDetectorDistance_mm(); + prepass_detector_geometry_.reset(); + prepass_rotation_scale_.reset(); + prepass_result_.reset(); + force_rotation_result_.reset(); + bragg_adaptive_.reset(); + beam_center_same_lattice_alternative_.reset(); + + logger.Info("Rotation two-pass geometry post-refinement: first pass (header geometry) -> {}_01_*", + base_prefix); + // Keep an empty prefix empty (the "compute stats, persist nothing" mode): the "_01" suffix would + // make it non-empty and RunPipeline would then write stray _01_* files despite no output wanted. + config_.output_prefix = base_prefix.empty() ? base_prefix : base_prefix + "_01"; + auto pass1 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true); + pass1.pass_number = 1; + pass1.pass_count = 2; + if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } + + // Beam centre as a LATTICE hypothesis, decided on the merge. The first pass indexed both + // centres, found the same cell with a different metric symmetry read off it, and measured that + // the centre disagreement is by itself large enough to have decided that symmetry (see the + // beam-centre check in RunPipeline). The class is handed to GeometryRefiner, which holds the + // cell to it for the rest of the run, so the wrong one costs an axis - refined, predicted, + // integrated and merged with an equality the crystal does not have. + // + // The disagreement is two-sided, so there is no centre to prefer: over the corpus runs that + // reach here the file's centre finds the higher class about as often as the measured one does. + // Nothing before integration separates them either - both index the same frames, and the frame + // count is the one statistic the check already refuses to arbitrate on. So run the first pass + // again at the other centre and keep whichever arm MERGES better, each judged on its own search + // merge: P1, whole range, before any correction surface, which is the merge the two passes are + // already compared in at the end of this function. Because the alternative arm is a whole first + // pass it gets the short-axis pass and the geometry post-refinement of its own, which the + // question needs - on a crystal with an axis under the FFT floor the check runs before that + // pass and compares supercells. + // + // Adopting a different metric moves the run off the geometry it was given, so it has to be + // decisively better rather than merely ahead. + // + // The same two arms also answer a second question the first pass cannot: where neither centre + // indexes the whole spot list and both index the same lattice once the pass reads less of each + // frame, which centre is right (see the beam-centre check). The measured arm then starts from + // the depth the check found it at. + // + // And a third: where both centres index the same lattice in the same class but lie further + // apart than the geometry absorbs, the centre can be off by about one reflection spacing - the + // same cell, every reflection one index out over part of the sweep - and only the merge shows + // it. Running the first pass twice on every such run would cost a pass on half of them, so it + // is run where the file's centre merged inconsistently: the lowest-resolution shell of the + // search merge, where every reflection is strong and noise cannot pull the correlation down, + // reads its half-sets as disagreeing. Measured on a pseudo-tetragonal crystal whose two centres + // lie 10.5 px apart, one c* spacing being 11.3 px: 0.610 at the file's centre and 0.997 at the + // measured one, R_meas of the final merge 0.58 against 0.13. + constexpr double MIN_INNER_CC_HALF = 0.9; + if (!cancelled_ && !beam_center_alternative_ && beam_center_same_lattice_alternative_ + && pass1.search_merge_inner_cc_half < MIN_INNER_CC_HALF) { + logger.Warning("Beam centre check: the file's centre merges inconsistently - CC1/2 {:.3f} in " + "the lowest-resolution shell of the search merge, where every reflection is " + "strong - and the measured centre, which indexes the same lattice, is further " + "away than this geometry absorbs, so the run tries it", + pass1.search_merge_inner_cc_half); + beam_center_alternative_ = beam_center_same_lattice_alternative_; + } else if (beam_center_same_lattice_alternative_) + logger.Info("Beam centre check: the file's centre merges consistently - CC1/2 {:.3f} in the " + "lowest-resolution shell of the search merge - so the measured centre, which " + "indexes the same lattice, is not tried", pass1.search_merge_inner_cc_half); + beam_center_same_lattice_alternative_.reset(); + if (!cancelled_ && beam_center_alternative_) { + const auto alt_center = *beam_center_alternative_; + beam_center_alternative_.reset(); + logger.Info("Beam centre check: running the first pass again at the measured centre " + "({:.2f},{:.2f}) and judging the two on their merges", alt_center[0], + alt_center[1]); + note(fmt::format("the first pass runs again at the measured beam centre ({:.1f}, {:.1f})", + alt_center[0], alt_center[1])); + const DiffractionExperiment file_experiment = experiment_; + const auto file_spot_finding = config_.spot_finding; + const auto file_mosaicity = prepass_mosaicity_; + const auto file_geometry = prepass_detector_geometry_; + const auto file_result = prepass_result_; + + 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(); + + 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. + 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; + // 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 + // FFT floor the check compared two supercells the short-axis pass then threw away. Where + // both arms now hold the same lattice in the same class, the centre did not decide the + // symmetry after all: the two merges are two integrations of ONE hypothesis, their CC1/2 + // difference is the noise between them, and it is not asked. Same lattice means what it + // means in the check: same class and primitive volumes within 2 %. + // + // What is left is which of the two GEOMETRIES is the more accurate, answered on the signal + // each arm measured (MeasuredCentreWins). The post-refinement cannot be left to settle it: + // it moves the beam by what the spots ask of a lattice already fitted at that centre, which + // on the crystal that showed it was 0.15 px of a 12.6 px error. + const bool arms_agree = ArmsHoldSameLattice(pass1, alt); + if (arms_agree) + logger.Info("Beam centre check: after their whole first passes both centres hold the " + "same lattice - {}-centred {} - so the centre did not decide the metric " + "symmetry; the two are judged on the signal each measured, {} reflections at " + "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); + 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", + pass1.rotation_lattice_type->centering, + 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); + if (!measured_wins.empty()) { + const auto cell_text = [](const ProcessResult &r) { + if (!r.consensus_cell) + return std::string("no cell"); + const auto &c = *r.consensus_cell; + return fmt::format("{:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}", + c.a, c.b, c.c, c.alpha, c.beta, c.gamma); + }; + logger.Warning("Beam centre check: {}, on {} against {} - so the run adopts the measured " + "centre and the lattice it finds", measured_wins, cell_text(alt), + cell_text(pass1)); + 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); + 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; + } + // 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. + beam_center_alternative_.reset(); + beam_center_alternative_state_.reset(); + beam_center_same_lattice_alternative_.reset(); + if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } + } + + // Undo the pre-pass goniometer shift (a sub-range or strided run moves start/increment inside + // the pass), but KEEP an axis direction the pass changed. The rotation-axis sign rescue runs + // inside the pass and adopts the opposite axis when the file's own indexes almost nothing; + // restoring the snapshot wholesale threw that away, handing pass 2 the sign that failed and + // then applying pass 1's fitted rotation scale to the axis it was not fitted on. It + // self-healed - pass 2's rescue fires again and a wrong sign indexes almost no frames - but + // only by repeating a whole first pass, and only while "a wrong sign always scores below + // half" holds, which is an invariant nothing states or tests. + if (gonio_snapshot) { + GoniometerAxis restored = *gonio_snapshot; + if (const auto after_pass1 = experiment_.GetGoniometer()) + restored.Axis(after_pass1->GetAxis()); + experiment_.Goniometer(restored); + } + // A measured spot budget is the pass's own, not the run's: give the second pass the same list to + // measure from, so the two passes cannot ratchet each other down. The canonical pass then reads + // the budget off the refined geometry, where the spots that do lie on the lattice actually do. + experiment_.MaxSpotCount(max_spot_count_snapshot); + + // Apply the post-refined detector geometry for the second pass, keeping the header geometry so + // the run can go back to it if the refined pass turns out worse (see the quality guard below). + const std::array header_geometry{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), + experiment_.GetDetectorDistance_mm(), + experiment_.GetPoniRot1_rad(), + experiment_.GetPoniRot2_rad()}; + if (prepass_detector_geometry_) { + const auto &g = *prepass_detector_geometry_; + experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) + .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); + } + // The probe passes below and the lattice arms after them exist ONLY to measure - see there - so + // the count of passes the run made has to include them even though they write nothing. + int arm_passes = 0; + + // The goniometer rotation scale: the angles in the file are the commanded ones, so a stage that + // turned at the wrong rate is a geometry error like any other - but one pass 1's fit reads only + // in part (see WalkRotationScale), so it is walked to the fit's fixed point here and adopted + // only where the validation frames of the whole sweep prefer it. The walk runs at the detector + // geometry the second pass will, so every hypothesis, the stored angles included, is scored + // there, and at that geometry only. Its passes + // only measure - an indexing probe stops once the lattice is scored, a refit stops once the + // post-refinement has measured - and they leave nothing behind: the experiment, the geometry + // the second pass is to run at and pass 1's mosaicity (a width in degrees fitted against the + // stored angles) are put back when the walk ends. + std::string rotation_scale_walk; + if (!cancelled_ && gonio_snapshot && pass1.post_refine && !config_.rotation_scale) { + step("Rotation-scale check"); + const DiffractionExperiment before_walk = experiment_; + const auto geometry_before_walk = prepass_detector_geometry_; + const auto mosaicity_before_walk = prepass_mosaicity_; + // A probe pass at scale k, run on r - this run, or a copy of it (see below). + const auto probe_on = [&](Rugnux &r, RugnuxObserver *obs, float k, bool index_only) { + r.experiment_ = before_walk; + r.experiment_.Goniometer(ScaleRotation(*before_walk.GetGoniometer(), k)); + r.prepass_rotation_scale_ = k; + r.prepass_mosaicity_.clear(); + r.indexing_probe_only_ = index_only; + r.postrefine_probe_ = r.postrefine_probe_only_ = !index_only; + ProcessResult res; + try { + res = r.RunPipeline(obs, /*write_output=*/false, /*geometry_prepass=*/false); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + // Angles under which nothing indexes have scored nothing, which is what an empty + // result reads as. + logger.Info("Rotation scale {:.5f}: the probe pass did not complete ({})", k, e.what()); + } + r.indexing_probe_only_ = r.postrefine_probe_ = r.postrefine_probe_only_ = false; + return res; + }; + const auto probe = [&](float k, bool index_only) { + ++arm_passes; + return probe_on(*this, observer, k, index_only); + }; + // The walk always indexes at the stored angles first and at the pass-1 fit next (see + // WalkRotationScale). Both are known now and each is its own pass, so the one at the fit + // runs at the same time, on a copy of the run as it stands before either. What a probe + // leaves behind for the next pass - the spot-finding settings a first-pass ladder rung + // adopts; the experiment is put back by the probe itself - is checked when the fit's turn + // comes: if the probe at the stored angles changed it, the copy started from a state the + // sequential walk would not have, and the probe is run again here. + const float first_fit = static_cast(pass1.post_refine->rotation_scale); + std::shared_ptr fit_copy; + std::future at_fit; + const SpotFindingSettings walk_start_spot_finding = config_.spot_finding; + if (first_fit != 1.0f && !cancelled_) { + fit_copy = std::make_shared(*this); + at_fit = std::async(std::launch::async, [&probe_on, fit_copy, first_fit] { + return probe_on(*fit_copy, nullptr, first_fit, true); + }); + } + const auto index_at = [&](float k) { + if (at_fit.valid() && k == first_fit) { + ProcessResult r = at_fit.get(); + const bool same_start = + config_.spot_finding.high_resolution_limit == walk_start_spot_finding.high_resolution_limit + && config_.spot_finding.measured_ring_q_recipA == walk_start_spot_finding.measured_ring_q_recipA; + if (same_start) { + ++arm_passes; + config_.spot_finding = fit_copy->config_.spot_finding; + for (const auto &m : fit_copy->first_pass_memo_) + AddFirstPassMemo(m); + return r.validation_evidence; + } + } + return probe(k, true).validation_evidence; + }; + const auto refit_at = [&](float k) -> std::optional { + const ProcessResult r = probe(k, false); + if (!r.post_refine) + return std::nullopt; + logger.Info("Rotation scale {:.5f}: the post-refinement there fits {:.5f} on top of it, " + "held-out residual {:.3e}", k, r.post_refine->rotation_scale, + r.post_refine->held_out_before); + return r.post_refine->rotation_scale; + }; + constexpr int MAX_ROTATION_SCALE_ROUNDS = 8; + const RotationScaleWalk walk = WalkRotationScale(pass1.post_refine->rotation_scale, index_at, + refit_at, MAX_ROTATION_SCALE_ROUNDS); + experiment_ = before_walk; + prepass_detector_geometry_ = geometry_before_walk; + prepass_mosaicity_ = mosaicity_before_walk; + prepass_rotation_scale_.reset(); + if (!walk.trail.empty()) { + rotation_scale_walk = fmt::format( + "goniometer rotation scale walked from the pass-1 fit {:.5f}, validation spots on " + "the lattice beyond chance by scale: {} - {}", pass1.post_refine->rotation_scale, + walk.trail, walk.scale == 1.0f ? "the stored angles stand" + : fmt::format("{:.5f} adopted", walk.scale)); + logger.Info("Two-pass: {}", rotation_scale_walk); + } + note(walk.scale == 1.0f ? std::string("the stored rotation angles stand") + : fmt::format("rotation scale {:.5f} adopted", walk.scale)); + if (walk.scale != 1.0f) { + prepass_rotation_scale_ = walk.scale; + pass1.post_refine->rotation_scale = walk.scale; + pass1.post_refine->rotation_scale_suspect = true; + logger.Warning("Goniometer rotation scale {:.5f} ({:+.2f} %): the stage turned {} than the " + "angles stored in the file, which are the COMMANDED values. The second pass " + "integrates at the corrected angles; the fault is in the hardware and should " + "be fixed there", walk.scale, 100.0 * (walk.scale - 1.0), + walk.scale > 1.0f ? "further" : "less far"); + } + } + + // ... and apply the adopted scale, on the same measure-then-re-integrate footing as the geometry. + if (prepass_rotation_scale_ && gonio_snapshot) { + experiment_.Goniometer(ScaleRotation(*experiment_.GetGoniometer(), *prepass_rotation_scale_)); + // The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has + // just stopped using, and the override can only ever raise the second pass's own estimate (it + // takes the larger of the two). Carrying it over would hold the second pass at the rocking + // width the uncorrected angles produced - the correction half-applied. Drop it and let the + // second pass fit its own. + prepass_mosaicity_.clear(); + logger.Info("Two-pass: goniometer rotation scale {:.5f} will drive the second integration pass " + "(oscillation {:.4f} -> {:.4f} deg per image)", *prepass_rotation_scale_, + gonio_snapshot->GetIncrement_deg(), + gonio_snapshot->GetIncrement_deg() * *prepass_rotation_scale_); + } + // The integration radius the pre-scan measured from the crystal's own spot width, applied now + // that the geometry is fitted. The pre-pass ran at the radius the run started with (see + // bragg_adaptive_), so everything pass 2 consumes from it - the post-refined distance and beam, + // the space group, the smoothed mosaicity - is what the fixed radius produces, and the widened + // disk changes only the intensities of the canonical pass. + if (bragg_adaptive_) { + experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_); + logger.Info("Two-pass: the canonical pass integrates at the measured spot width, " + "r1={:.1f} r2={:.1f} r3={:.2f}", bragg_adaptive_->GetR1(), + bragg_adaptive_->GetR2(), bragg_adaptive_->GetR3()); + } + + // Let every refined pass from here on ask the quality guard's question of itself, so one that + // is going to be thrown away can skip the report and the files it would have superseded. + quality_guard_pass1_ = pass1; + config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix) + // The canonical pass measures the post-refinement again where the pre-pass moved the geometry + // (it applies nothing), so that the walk below has something to decide on. + postrefine_probe_ = prepass_detector_geometry_.has_value(); + + const auto current_geometry = [&] { + return std::array{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), + experiment_.GetDetectorDistance_mm(), + experiment_.GetPoniRot1_rad(), experiment_.GetPoniRot2_rad()}; + }; + const auto set_geometry = [&](const std::array &g) { + experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) + .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); + }; + + // The metric symmetry, asked of the spots rather than of the tolerance that admitted it. The + // Bravais class is chosen by a walk that reads two axes as EQUAL when they agree to a fixed + // relative tolerance, and the class carrying that equality is then imposed everywhere below: + // the cell is refined with a = b, reflections are predicted from it and integrated at the + // predicted positions. On a small cell the tolerance is far wider than the positions resolve - + // a 5.8 A axis pair 2.4 % apart is read as tetragonal, and the mean the constraint substitutes + // mispredicts the outermost reflections by a good ten pixels against an integration disc of + // four - so a genuinely orthorhombic crystal is integrated as tetragonal and every decision + // downstream, up to the space group, is read off the wreckage. + // + // So the equality is measured rather than assumed. The indexer already refines the same + // candidate a second time with nothing held, and pass 1 reports how far that free refinement + // leaves the two axes apart (lattice_length_split). A relative split of eps displaces a + // reflection at radius r by eps/2 * r pixels, so the question "do these data resolve it" has + // an answer in the units the integration works in: where that displacement at the far corner + // of the detector stays inside the integration disc, imposing the equality moves nothing out + // of its own box and the higher symmetry is kept with no pass run. Where it does not, both + // hypotheses are run - the promoted class, and the class the same walk carries when it is + // granted no length equality at all - and whichever REALISES the lower held-out residual on a + // re-integrated, re-indexed pass is kept, the + // promoted class on a tie, because it is the one that asserts more. Both arms are probe-only + // passes (postrefine_probe_only_): the residual is all the comparison reads and it is in hand + // before the scaling engine is built, so neither arm pays for a merge, a space-group search or + // a set of correction surfaces, and only the class that wins is ever run as a canonical pass. + std::string lattice_arm; // the arm decision, for the pass decision text + // + // An ANGLE promotion is the same question: the walk holds a monoclinic or orthorhombic angle at + // 90 when the unrefined cell puts it within LATTICE_SEARCH_ANGLE_TOLERANCE_DEG of it, and an + // angle error of delta displaces a reflection at radius r by up to delta/2 * r pixels. Its + // demoted arm is the indexer's free refinement itself; the space group, and with it any + // lower class the data do carry, is then found from the intensities as for any triclinic + // indexing. Measured: a monoclinic crystal whose beta sits 1.4 deg from 90 indexed every + // validation frame as primitive orthorhombic, so the per-frame guard never fired. + const bool length_arm = pass1.lattice_length_split.has_value(); + const bool angle_arm = !length_arm && pass1.lattice_angle_departure_deg.has_value(); + if (!cancelled_ && (length_arm || angle_arm) && pass1.rotation_lattice_type) { + const double split = length_arm ? *pass1.lattice_length_split + : *pass1.lattice_angle_departure_deg * PI / 180.0; + const double dx = std::max(experiment_.GetBeamX_pxl(), + static_cast(experiment_.GetXPixelsNum()) - experiment_.GetBeamX_pxl()); + const double dy = std::max(experiment_.GetBeamY_pxl(), + static_cast(experiment_.GetYPixelsNum()) - experiment_.GetBeamY_pxl()); + const double r_corner_px = std::hypot(dx, dy); + const double miss_px = 0.5 * split * r_corner_px; + const double disc_px = experiment_.GetBraggIntegrationSettings().GetR1(); + if (miss_px > disc_px) { + if (length_arm) + logger.Info("Metric symmetry: {} was promoted on a = b, but the free refinement leaves " + "the two axes {:.2f} % apart - {:.1f} px of misprediction at the far corner " + "against an integration disc of {:.1f} px. Both classes are re-integrated " + "and judged on the held-out residual each realises", + gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), + 100.0 * split, miss_px, disc_px); + else + logger.Info("Metric symmetry: {} holds its angles at 90 deg, but the free refinement " + "leaves one {:.2f} deg from it - {:.1f} px of misprediction at the far " + "corner against an integration disc of {:.1f} px. Both classes are " + "re-integrated and judged on the held-out residual each realises", + gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), + *pass1.lattice_angle_departure_deg, miss_px, disc_px); + step("Metric-symmetry check"); + note(fmt::format("{} re-integrated with and without its {} equality", + gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), + length_arm ? "a = b" : "90 deg angle")); + // The arms are compared on the post-refinement's held-out residual, so they have to + // measure one: a canonical pass only does where the pre-pass moved the geometry. + const bool probe_was = postrefine_probe_; + postrefine_probe_ = true; + postrefine_probe_only_ = true; + const auto promoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + ++arm_passes; + IndexingSettings demoted_settings = experiment_.GetIndexingSettings(); + demoted_settings.LatticeClassDemoted(length_arm); + experiment_.ImportIndexingSettings(demoted_settings); + lattice_angles_demoted_ = angle_arm; + const auto demoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + ++arm_passes; + postrefine_probe_only_ = false; + postrefine_probe_ = probe_was; + const double realised_promoted = promoted_arm.post_refine ? promoted_arm.post_refine->held_out_before : NAN; + const double realised_demoted = demoted_arm.post_refine ? demoted_arm.post_refine->held_out_before : NAN; + const auto class_text = [](const ProcessResult &r) { + return r.rotation_lattice_type + ? fmt::format("{}-centred {}", r.rotation_lattice_type->centering, + gemmi::crystal_system_str(r.rotation_lattice_type->crystal_system)) + : std::string("no lattice"); + }; + if (!cancelled_ && promoted_arm.post_refine && demoted_arm.post_refine + && HeldOutResidualFell(*promoted_arm.post_refine, *demoted_arm.post_refine)) { + logger.Warning("Metric symmetry: re-integrated as {} the held-out residual is {:.3e} " + "against {:.3e} as {} - the equality is not this crystal's, and the " + "run goes on without it", class_text(demoted_arm), realised_demoted, + realised_promoted, class_text(promoted_arm)); + note(fmt::format("the equality dropped: {} fits better", class_text(demoted_arm))); + lattice_arm = fmt::format( + "the metric symmetry was promoted to {} on {} the free refinement leaves " + "{:.2f} {}; re-integrated as {} the held-out residual is {:.3e} against " + "{:.3e}, so the equality was dropped", class_text(promoted_arm), + length_arm ? "axes" : "an angle", length_arm ? 100.0 * split : split * 180.0 / PI, + length_arm ? "% apart" : "deg from 90", + class_text(demoted_arm), realised_demoted, realised_promoted); + } else if (!cancelled_) { + note(fmt::format("the metric symmetry stands: {}", class_text(promoted_arm))); + IndexingSettings promoted_settings = experiment_.GetIndexingSettings(); + promoted_settings.LatticeClassDemoted(false); + experiment_.ImportIndexingSettings(promoted_settings); + lattice_angles_demoted_ = false; + // An arm that indexed nothing, or post-refined nothing, measures no residual. The + // metric symmetry stands either way - that is the arm whose absence of evidence is + // not evidence against it - but the run must not report a comparison it could not + // make: written as a number, a missing residual came out as "realised a held-out + // residual of nan against 1.234e-03, so the metric symmetry stands". + const bool compared = std::isfinite(realised_promoted) && std::isfinite(realised_demoted); + if (compared) { + logger.Info("Metric symmetry: re-integrated as {} the held-out residual is {:.3e}, " + "which {} does not beat ({:.3e}) - the higher symmetry stands", + class_text(promoted_arm), realised_promoted, class_text(demoted_arm), + realised_demoted); + lattice_arm = fmt::format( + "{} was run as well and realised a held-out residual of {:.3e} against {:.3e} " + "as {}, so the metric symmetry stands", class_text(demoted_arm), + realised_demoted, realised_promoted, class_text(promoted_arm)); + } else { + const char *which = !std::isfinite(realised_demoted) && !std::isfinite(realised_promoted) + ? "neither arm measured a held-out residual" + : !std::isfinite(realised_demoted) + ? "the demoted arm measured no held-out residual" + : "the promoted arm measured no held-out residual"; + logger.Info("Metric symmetry: {}, so the two could not be compared - the higher " + "symmetry stands", which); + lattice_arm = fmt::format( + "{} was run as well, but {}, so the two could not be compared and the " + "metric symmetry stands untested", class_text(demoted_arm), which); + } + } + } else { + logger.Info("Metric symmetry: {} - the free refinement departs from it by {:.2f} {}, " + "{:.1f} px at the far corner, inside the integration disc of {:.1f} px, so " + "these data do not resolve the difference", + gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), + length_arm ? 100.0 * split : split * 180.0 / PI, length_arm ? "%" : "deg", + miss_px, disc_px); + } + } + + logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*", + base_prefix); + speculate_geometry_probe_ = true; // a walk that follows starts with an indexing probe + rerun_if_starved_ = true; + auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + rerun_if_starved_ = false; + speculate_geometry_probe_ = false; + + // Walk the geometry to the fit's fixed point. The post-refinement is fitted to reflections + // integrated at the very geometry it is correcting, so one fit reaches only as far as those + // reflections see: on a header whose distance is several percent long, the pre-pass absorbs the + // error into an inflated cell, and from that self-consistent pair the fit finds a fraction of it + // and stops. The canonical pass has just re-integrated and re-indexed at what it committed, so + // asking the same fit there is the next step of the same walk - and that re-indexing is also the + // RATIFICATION the old size bound stood in for: a second lattice does not index better at a + // geometry it did not come from, a wrong header does, and a pass that lost the crystal is thrown + // out by the quality guard below, which returns the run to the header geometry. + // + // What starts the walk outright is a move of more than one step: that is the fit saying the + // geometry is somewhere else entirely rather than a fraction of a percent away. A smaller move + // has to be preferred by the validation frames first (below), so a run whose fit settles beside + // the geometry in hand costs two indexing probes and no pass. Once the run IS walking + // it keeps walking because the steps do not get smaller in proportion to what is left: measured + // on a header 3.9 % long, the first round took a third of the error and each of the next took a + // third of the rest. + // + // But a round is kept only for what it REALISES, not for what the fit predicts: the held-out + // residual of the pass re-integrated and re-indexed at the new geometry has to fall below the + // best round's by more than the cross-validation noise of the two (HeldOutResidualFell). With + // the walk allowed twelve rounds, every round's fit predicted a gain the re-integration did not + // show - the realised residual stayed flat, or passed a minimum and rose again, while the + // distance crept on along the distance/cell-scale degeneracy - so "the fit still commits" never + // ended a walk, and the cap decided where on the drift the run stopped. A round that does not + // pay ends the walk and the run goes back to the best round it reached. So does a fit whose + // committed cell the re-indexing pushes back the same way twice running (ReindexPushesCellBack): + // the fit keeps asking for a cell the indexer at the fit's own geometry does not return, and + // nothing ratifies it. A walk still paying when the rounds run out keeps its best round - the + // last - and says it did not converge; the quality guard below judges it against the header + // pass as it judges every refined pass. + // + // The residual is not the only thing a round realises. It is the mean over every reflection, + // measured by centroids taken inside a disc centred on the prediction - so they follow the + // prediction part of the way - and dominated by the low-resolution reflections, where a distance + // and the cell scale that compensates it move every spot alike. What does see such an error is + // the re-indexing: the share of the validation spots that land on the lattice, which the + // high-resolution spots lose first. Measured on a crystal whose canonical pass ran 0.6 % long, + // the round at the fit's distance put 70.4 % of the validation spots on the lattice against + // 58.1 %, while reading the residual only 0.74 of its noise lower - and the run the residual + // stopped kept the long distance and lost 0.05 A of resolution. So a round is kept when EITHER + // realises a gain beyond its noise: the residual (HeldOutResidualFell), or the validation + // frames (ValidationEvidencePrefers). + // + // And a move of less than one step no longer settles the geometry by being small: it is asked + // of the validation frames too, but cheaply first - an indexing probe at the fit's geometry + // against one at the geometry in hand, each scoring the lattice and stopping there - and only a + // move the probes prefer pays for a re-integrated round. A move of a step or more is run as a + // round directly, as it always was, and so is every round of the walk it starts. + const auto on_lattice_pct = [](const ValidationSpotEvidence &e) { + return 100.0 * static_cast(e.on_lattice - e.by_chance) + / static_cast(std::max(1, e.spots)); + }; + // An indexing probe at a geometry: the validation evidence of the lattice indexed there, and + // nothing else - the pass stops once the lattice is scored, and the experiment is put back. + const auto index_at_geometry = [&](const std::array &g) { + const DiffractionExperiment before_probe = experiment_; + const bool probe_was = postrefine_probe_; + set_geometry(g); + CollectSpeculativeProbe(); + indexing_probe_only_ = true; + postrefine_probe_ = false; + ProcessResult r; + try { + r = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Info("Two-pass: the indexing probe at distance {:.3f} mm did not complete ({})", + g[2], e.what()); + } + indexing_probe_only_ = false; + postrefine_probe_ = probe_was; + experiment_ = before_probe; + ++arm_passes; + return r.validation_evidence; + }; + constexpr int MAX_GEOMETRY_ROUNDS = 8; + int geometry_rounds = 0; // passes run at a geometry the walk moved to + int walk_passes = 0; // every pass the walk ran, the return to its best round included + int best_round = 0; + std::array best_geometry = current_geometry(); + PostRefineResult best_fit; // the fit measured in the best round's pass: its realised residual + int last_pushback = 0; + bool geometry_not_converged = false; + std::string walk_stop; // why the walk ended, where anything but the fit settling ended it + std::string walk_rounds; // distance / realised residual per round, for the pass decision + if (pass2.post_refine) { + best_fit = *pass2.post_refine; + walk_rounds = fmt::format("{:.3f} mm / {:.3e} / {:.1f}%", best_geometry[2], best_fit.held_out_before, + on_lattice_pct(pass2.validation_evidence)); + } + ValidationSpotEvidence best_evidence = pass2.validation_evidence; + std::optional probe_at_best; // the indexing probe at best_geometry, once run + bool large_walk = false; // a round of this walk was started by a move of a step or more + // What the last round started from: a pass can leave the run changed behind it (a spot budget + // it shortened, an axis sign it flipped), so a round is written + // from here, not from where its own probe left the run. + DiffractionExperiment last_round_start = experiment_; + bool walk_shown = false; + while (!cancelled_ && prepass_detector_geometry_ + && pass2.post_refine && pass2.post_refine->detector_refined) { + if (!walk_shown) { + step("Geometry check"); + walk_shown = true; + } + if (geometry_rounds == MAX_GEOMETRY_ROUNDS) { + geometry_not_converged = true; + walk_stop = fmt::format("the fit still moved the geometry after {} rounds, every one of " + "which lowered the realised residual", MAX_GEOMETRY_ROUNDS); + logger.Warning("Two-pass: {} - the walk did not converge; the run keeps its last round " + "and the quality guard judges it against the header pass", walk_stop); + break; + } + const PostRefineResult fit = *pass2.post_refine; + const std::array g = *prepass_detector_geometry_; + large_walk = large_walk || fit.large_move; + if (!large_walk) { + if (!probe_at_best) + probe_at_best = index_at_geometry(best_geometry); + const ValidationSpotEvidence probe_at_fit = index_at_geometry(g); + const bool prefers = ValidationEvidencePrefers(*probe_at_best, probe_at_fit); + logger.Info("Two-pass: the post-refinement moves the geometry by less than a step (distance " + "{:.3f} -> {:.3f} mm, beam {:.2f},{:.2f} -> {:.2f},{:.2f} px); indexing probes put " + "{:.1f}% of the validation spots on the lattice beyond chance there against {:.1f}% " + "at the geometry in hand - {}", fit.distance_before_mm, fit.distance_after_mm, + fit.beam_x_before_px, fit.beam_y_before_px, fit.beam_x_after_px, fit.beam_y_after_px, + on_lattice_pct(probe_at_fit), on_lattice_pct(*probe_at_best), + prefers ? "worth a round" : "the geometry in hand stands"); + if (!prefers) + break; + probe_at_best = probe_at_fit; + } + logger.Info("Two-pass: the post-refinement at the adopted geometry still moves it " + "(distance {:.3f} -> {:.3f} mm) - re-integrating and re-indexing at what it " + "asks for (round {} of at most {})", fit.distance_before_mm, + fit.distance_after_mm, geometry_rounds + 1, MAX_GEOMETRY_ROUNDS); + set_geometry(g); + last_round_start = experiment_; + // A round is a probe pass: the walk reads only its fit and its validation evidence, both in + // hand before the merge, so the files are written once, below, at the round the walk keeps. + speculate_geometry_probe_ = !large_walk; + postrefine_probe_only_ = true; + pass2 = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + postrefine_probe_only_ = false; + speculate_geometry_probe_ = false; + ++geometry_rounds; + ++walk_passes; + const double realised = pass2.post_refine ? pass2.post_refine->held_out_before : NAN; + walk_rounds += fmt::format(", {:.3f} mm / {:.3e} / {:.1f}%", g[2], realised, + on_lattice_pct(pass2.validation_evidence)); + const bool residual_fell = pass2.post_refine && HeldOutResidualFell(best_fit, *pass2.post_refine); + const bool evidence_prefers = ValidationEvidencePrefers(best_evidence, pass2.validation_evidence); + if (!pass2.post_refine || !(residual_fell || evidence_prefers)) { + walk_stop = fmt::format( + "round {} lowered neither the realised held-out residual by more than its noise " + "({:.3e} against {:.3e}, standard errors {:.1e} and {:.1e}) nor the share of " + "validation spots off the lattice ({:.1f}% on it beyond chance against {:.1f}%)", + geometry_rounds, realised, best_fit.held_out_before, + pass2.post_refine ? pass2.post_refine->held_out_before_se : NAN, + best_fit.held_out_before_se, on_lattice_pct(pass2.validation_evidence), + on_lattice_pct(best_evidence)); + logger.Info("Two-pass: {}", walk_stop); + break; + } + best_round = geometry_rounds; + best_geometry = g; + best_fit = *pass2.post_refine; + best_evidence = pass2.validation_evidence; + const int pushback = ReindexPushesCellBack(fit, *pass2.post_refine); + if (pushback != 0 && pushback == last_pushback) { + walk_stop = fmt::format( + "re-indexing pushed the cell length {} back against the fit for the second " + "round running, so nothing ratifies the cell the fit asks for", + "abc"[std::abs(pushback) - 1]); + logger.Info("Two-pass: {} - the walk stops at round {}", walk_stop, geometry_rounds); + break; + } + last_pushback = pushback; + } + postrefine_probe_ = false; + if (walk_shown) + note(geometry_rounds == 0 ? std::string("the geometry in hand stands") + : fmt::format("geometry walked {} round{}, round {} kept: distance {:.2f} mm{}", geometry_rounds, + geometry_rounds == 1 ? "" : "s", best_round, best_geometry[2], + best_round > 0 ? " - the files are written again there" : "")); + // The rounds wrote nothing, so the round the walk kept is run again to write the canonical files. + if (!cancelled_ && geometry_rounds > 0) { + if (best_round == geometry_rounds) { + logger.Info("Two-pass: writing round {} of the geometry walk (distance {:.3f} mm)", + best_round, best_geometry[2]); + experiment_ = last_round_start; + } else { + logger.Info("Two-pass: going back to round {} of the geometry walk (distance {:.3f} mm), " + "the last one that lowered the realised residual", best_round, best_geometry[2]); + set_geometry(best_geometry); + } + rerun_if_starved_ = true; + pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + rerun_if_starved_ = false; + ++walk_passes; + } + + // The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever + // result is returned - it is a result of the run, not of the pass that happened to fit it. Where + // the walk took rounds they are ONE geometry decision, so what the report names is where the + // geometry started and where the run left it, not the first step of the way there. + if (best_round > 0 && pass1.post_refine) { + pass1.post_refine->beam_x_after_px = best_geometry[0]; + pass1.post_refine->beam_y_after_px = best_geometry[1]; + pass1.post_refine->distance_after_mm = best_geometry[2]; + // A first pass that held the header distance and a walk that then left it: the run's + // decision is the walk's. + pass1.post_refine->distance_held = false; + } + pass2.post_refine = pass1.post_refine; + pass2.pass_number = 2 + arm_passes + walk_passes; + pass2.pass_count = 2 + arm_passes + walk_passes; + if (geometry_rounds > 0) + pass2.pass_decision = fmt::format( + "post-refined geometry walked {} round{} and kept round {}, detector distance {:.3f} mm " + "against the header's {:.3f}: {}; distance / realised held-out residual / validation spots on the lattice beyond chance by round: {}", + geometry_rounds, geometry_rounds == 1 ? "" : "s", best_round, best_geometry[2], + pass1.post_refine ? pass1.post_refine->distance_before_mm : 0.0, + walk_stop.empty() ? "the fit at the kept geometry commits no further move" : walk_stop, + walk_rounds); + + // r2 = r1 + 2 is also the inner edge of the background ring, so widening the signal disk pushes + // that ring out into the neighbours; a reflection whose ring is left with five or fewer clean + // pixels has no background and is dropped whole. On a pattern dense enough in three dimensions + // at once that costs a quarter of the observations - and NOTHING available before integrating + // predicts it: the crystal it happens to has a mid-table predicted spot spacing, wider than + // crystals that survive r1 = 12. What does measure it is the integrator counting its own drops. + // The pass that integrates at the widened radius is now the canonical one, so the measurement + // and the remedy are both here: re-run it at the radius that keeps the data, the same remedy the + // two guards below use and, like them, a pass only on a crystal that was going to be wrong. + // Asked of every pass that can end the run - the two guards below re-run at whatever radius is + // in hand, and a starved refined pass has written nothing (see the superseded flag in + // RunPipeline) - so each is followed by it. After the first re-run the radius is the fixed one + // and this never fires again. + const auto rerun_if_starved = [&](ProcessResult &pass, const std::string decision) { + if (cancelled_ || !bragg_adaptive_ || !bkg_starved_fraction_ + || *bkg_starved_fraction_ <= spot_width::BKG_STARVED_MAX_FRACTION) + return; + logger.Info("Two-pass: at the widened integration radius r1={:.1f} the neighbouring " + "reflections left {:.2f}% of the predicted reflections without a background " + "ring (bound {:.2f}%), so this pass is re-run at r1={:.1f} r2={:.1f} r3={:.2f}", + bragg_adaptive_->GetR1(), 100.0 * *bkg_starved_fraction_, + 100.0 * spot_width::BKG_STARVED_MAX_FRACTION, + bragg_before_adaptive_->GetR1(), bragg_before_adaptive_->GetR2(), + bragg_before_adaptive_->GetR3()); + const double starved = *bkg_starved_fraction_; + const float wide_r1 = bragg_adaptive_->GetR1(); + const int passes = pass.pass_count + 1; + experiment_.ImportBraggIntegrationSettings(*bragg_before_adaptive_); + bragg_adaptive_.reset(); + note(fmt::format("{:.1f}% of the background rings starved: the pass runs again at r1 {:.1f} px", + 100.0 * starved, bragg_before_adaptive_->GetR1())); + pass = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + pass.post_refine = pass1.post_refine; + pass.pass_number = passes; + pass.pass_count = passes; + pass.pass_decision = fmt::format( + "{}, re-run at the fixed integration radius: at r1={:.1f} neighbouring reflections " + "left {:.2f}% of the predicted reflections without a background ring", decision, wide_r1, + 100.0 * starved); + }; + rerun_if_starved(pass2, "post-refined geometry adopted"); + + // Supercell-collapse guard for the second pass. The de-novo lattice search can be BISTABLE on a + // periodically-modulated dataset: the tiny post-refine geometry move tips it out of the true-cell + // basin onto a spurious SUPERCELL (a large multiple of pass-1's cell volume), and the space group + // reused from pass-1 is then stamped onto that wrong-shaped cell - catastrophic. Pass 1 (header + // geometry, un-perturbed) found the correct lattice, so re-run the second pass with pass-1's FULL + // indexing result forced (its refined lattice + orientation + rotation axis + geometry - not just the + // cell, which loses the per-frame orientation and integrates worse): the re-run then matches de-novo + // quality but at the refined geometry. Only fires when pass 2 is markedly LARGER than pass 1 - a + // genuine de-novo demotion to a smaller primitive (the pseudo-symmetry case the second pass re-indexes + // de novo to catch) has pass2 <= pass1 and is left untouched. + if (!cancelled_ && prepass_result_ && pass1.consensus_cell && pass2.consensus_cell) { + // PRIMITIVE volumes, like the two guards in RunPipeline: a centred conventional cell is an + // exact integer multiple of its primitive one (C/I 2x, R 3x, F 4x), so two settings of the + // SAME lattice differ by that factor and comparing the conventional cells reads a mere + // change of setting as a supercell - which then forces pass 1's setting and, with it, its + // lower symmetry. Measured: every dataset this guard fired on was at an exact centring + // multiplicity, one of them an F-centred lattice held down to P1 at 3.99x. + const auto primitive_volume = [](const UnitCell &cell, std::optional centering) { + return std::abs(CrystalLattice(cell).ToPrimitive(centering.value_or('P')).CalcVolume()); + }; + const double v1 = primitive_volume(*pass1.consensus_cell, pass1.consensus_centering); + const double v2 = primitive_volume(*pass2.consensus_cell, pass2.consensus_centering); + if (v1 > 0.0 && v2 > 1.5 * v1) { + logger.Info("Two-pass: second-pass primitive cell volume {:.0f} A^3 is {:.2f}x pass-1's " + "({:.0f} A^3) - a spurious supercell; re-running the second pass with " + "pass-1's result forced", + v2, v2 / v1, v1); + note("a spurious supercell: the pass runs again on the pre-pass lattice"); + force_rotation_result_ = *prepass_result_; + // Pass 1's result carries pass 1's goniometer, and the per-image path takes its angles + // from there - forcing it whole would put the uncorrected angles back and silently undo + // the rotation scale for the re-run. + if (prepass_rotation_scale_ && force_rotation_result_->axis) + force_rotation_result_->axis = + ScaleRotation(*force_rotation_result_->axis, *prepass_rotation_scale_); + rerun_if_starved_ = true; + pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + rerun_if_starved_ = false; + pass2.post_refine = pass1.post_refine; + pass2.pass_number = 3; + pass2.pass_count = 3; + pass2.pass_decision = fmt::format( + "post-refined geometry adopted, re-run with pass-1's lattice forced: the de-novo " + "primitive cell volume {:.0f} A^3 was {:.2f}x pass-1's ({:.0f} A^3), a spurious " + "supercell", v2, v2 / v1, v1); + // Still forced: the re-run at the fixed radius is the same pass. + rerun_if_starved(pass2, pass2.pass_decision); + force_rotation_result_.reset(); + } + } + // The refined pass is normally the better answer, which is why it is the canonical output - + // but it is not guaranteed to be, and the guard that decides is RefinedPassIsWorse above. It + // has already been asked once, inside the pass itself, which is how a pass that is about to be + // thrown away knows not to write a report; this is where the answer is acted on. + if (!cancelled_ && pass1.has_merge_statistics && pass2.has_merge_statistics) { + logger.Info("Two-pass quality guard: the refined pass holds {} reflections at I/sigma >= 2 and " + "{} at <= -2 of {} in a {:.0f} A^3 primitive cell, the header-geometry pass {} and {} " + "of {} in {:.0f} A^3; low-order axial {} against {}", + pass2.search_merge_strong_reflections, pass2.search_merge_negative_reflections, + pass2.search_merge_reflections, PrimitiveCellVolume(pass2), + pass1.search_merge_strong_reflections, pass1.search_merge_negative_reflections, + pass1.search_merge_reflections, PrimitiveCellVolume(pass1), + pass2.search_merge_axial_reflections, pass1.search_merge_axial_reflections); + const std::string worse = RefinedPassIsWorse(pass1, pass2); + if (!worse.empty()) { + logger.Warning("Two-pass: the refined pass is worse than the header-geometry pass ({}) " + "- going back to the header geometry. The refined geometry did not help " + "this crystal.", worse); + // The refined pass has already written the canonical files, so re-run at the header + // geometry to replace them - the same remedy the supercell collapse above uses, and it + // only costs a pass on a crystal that was going to be wrong otherwise. + experiment_.BeamX_pxl(header_geometry[0]).BeamY_pxl(header_geometry[1]) + .DetectorDistance_mm(header_geometry[2]) + .PoniRot1_rad(header_geometry[3]).PoniRot2_rad(header_geometry[4]); + if (prepass_rotation_scale_ && gonio_snapshot) { + // The header rotation angles with it - but not the axis DIRECTION, which is not + // part of the geometry being reverted and may have been corrected by the sign + // rescue (see the restore before pass 2). + GoniometerAxis restored = *gonio_snapshot; + if (const auto current = experiment_.GetGoniometer()) + restored.Axis(current->GetAxis()); + experiment_.Goniometer(restored); + prepass_rotation_scale_.reset(); + pass1.post_refine->rotation_scale_suspect = false; + } + config_.output_prefix = base_prefix; + // This pass is the answer whatever it measures - the guard has had its one chance - + // so it must not conclude from its own search merge that it is going to be re-run. + quality_guard_pass1_.reset(); + // What the guard preferred is pass 1 - its lattice, orientation and axis at this very + // geometry - so that is what is re-integrated, not a fresh de-novo index. A de-novo + // re-index here is a new hypothesis nobody judged: measured, at the same centre where + // pass 1 indexed 56/60 validation frames it found 5/60, took the opposite axis sign and + // shipped a merge with R_meas in the hundreds of percent. + if (prepass_result_) + force_rotation_result_ = *prepass_result_; + rerun_if_starved_ = true; + note("the refined geometry is worse: the pass runs again at the header geometry"); + auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + rerun_if_starved_ = false; + redo.post_refine = pass1.post_refine; + redo.pass_number = pass2.pass_count + 1; + redo.pass_count = pass2.pass_count + 1; + redo.pass_decision = fmt::format( + "header geometry re-adopted: the post-refined pass was worse ({})", worse); + // Still forced: the re-run at the fixed radius is the same pass. + rerun_if_starved(redo, redo.pass_decision); + force_rotation_result_.reset(); + return redo; + } + if (pass2.pass_decision.empty()) + pass2.pass_decision = fmt::format( + "post-refined geometry adopted ({}, {})", SearchMergeCompletenessText(pass1, pass2), + SignalText(pass1, pass2)); + } + if (pass2.pass_decision.empty()) + pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_) + ? "post-refined geometry adopted" + : "the post-refinement committed no geometry change, so this pass reproduces the first"; + if (!lattice_arm.empty()) + pass2.pass_decision = lattice_arm + "; " + pass2.pass_decision; + if (!rotation_scale_walk.empty()) + pass2.pass_decision = rotation_scale_walk + "; " + pass2.pass_decision; + pass2.geometry_not_converged = geometry_not_converged; + return pass2; + } + return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); +}