diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index decacb39c..0ac770c0a 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1539,7 +1539,96 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { return result; } +namespace { + // How much better one pass's search merge has to be for the run to move off the other. Used by the + // beam-centre arbitration and by the cross-pass quality guard, both 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); + } + + // Pass 2 is normally the better answer, which is why it is the canonical output - but it is not + // guaranteed to be. Two ways it can be wrong: it merges more unique reflections than the cell it + // settled on can hold (completeness above 100% is arithmetically impossible and means the cell is + // wrong), or its CC1/2 collapses relative to pass 1. Both bounds are set where only a failure + // reaches them, so a normal run - where pass 2 is a little better - keeps pass 2. Returns why the + // refined pass is worse, in the words the run reports it in, or an empty string when it is not. + // + // 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 CC1/2 decides alone. 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. Two cheaper stand-ins were measured on that corpus and failed: + // pass 2's search-merge OBSERVATION count never drops below 90% of pass 1's even when pass 2 + // predicts 44% fewer partials, and pass 2's search-merge COMPLETENESS relative to pass 1's + // fires only on passes that reached FINER, because this merge is never resolution-cut and the + // possible list grows with the range - it reverted two healthy crystals and rescued none. + // + // 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 + // completeness and CC1/2 in different groups are not the same measurement. The search merge is + // P1 on both sides, over the full resolution range and before any correction surface is fitted + // - the one merge the two passes make in the same terms. + // + // 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 + // The third arm: the axial rows. CC1/2 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 CC1/2 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, which is exactly what a wedge-losing pass looks + // like from CC1/2 - better. 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 only decides where CC1/2 does not: a pass that wins on CC1/2 by more than the loss the arm + // above already tolerates has said something about the bulk of the data that a few dozen + // reflections cannot answer, and is kept. And it needs pass 1 to have measured enough of the + // rows for the comparison to mean anything - below eight there is no zone a screw could be + // claimed on either side, so the difference decides nothing. + constexpr double MIN_AXIAL_ROWS_KEPT = 0.75; + constexpr int64_t MIN_AXIAL_ROWS_TO_COMPARE = 8; + const double cc1 = pass1.search_merge_cc_half; + const double cc2 = pass2.search_merge_cc_half; + const int64_t axial1 = pass1.search_merge_axial_reflections; + const int64_t axial2 = pass2.search_merge_axial_reflections; + const bool lost_axial_rows = axial1 >= MIN_AXIAL_ROWS_TO_COMPARE + && axial2 < MIN_AXIAL_ROWS_KEPT * static_cast(axial1) + && cc2 < cc1 + MAX_CC_HALF_LOSS; + 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) + || cc2 < cc1 - MAX_CC_HALF_LOSS || lost_axial_rows; + if (!worse) + return {}; + return (cc2 < cc1 - MAX_CC_HALF_LOSS || !lost_axial_rows) + ? fmt::format("{}, CC1/2 before corrections {:.3f} vs {:.3f}", + SearchMergeCompletenessText(pass1, pass2), cc2, cc1) + : fmt::format("it holds {} of the {} low-order axial reflections the header-geometry " + "pass measured, and its CC1/2 before corrections {:.3f} vs {:.3f} does " + "not say it is better", axial2, axial1, cc2, cc1); + } +} + 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 @@ -1554,9 +1643,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { // 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_; - // How much better one pass's search merge has to be for the run to move off the other. Used - // twice here - by the beam-centre arbitration and by the pass-2 quality guard at the end. - constexpr double MAX_CC_HALF_LOSS = 0.05; const int64_t max_spot_count_snapshot = experiment_.GetDatasetSettings().GetMaxSpotCount(); prepass_detector_geometry_.reset(); prepass_rotation_scale_.reset(); @@ -1722,6 +1808,9 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { 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; logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*", base_prefix); config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix) @@ -1810,80 +1899,13 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { force_rotation_result_.reset(); } } - // Pass 2 is normally the better answer, which is why it is the canonical output - but it is not - // guaranteed to be. Two ways it can be wrong: it merges more unique reflections than the cell it - // settled on can hold (completeness above 100% is arithmetically impossible and means the cell is - // wrong), or its CC1/2 collapses relative to pass 1. Both bounds are set where only a failure - // reaches them, so a normal run - where pass 2 is a little better - keeps pass 2. - // - // 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 CC1/2 decides alone. 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. Two cheaper stand-ins were measured on that corpus and failed: - // pass 2's search-merge OBSERVATION count never drops below 90% of pass 1's even when pass 2 - // predicts 44% fewer partials, and pass 2's search-merge COMPLETENESS relative to pass 1's - // fires only on passes that reached FINER, because this merge is never resolution-cut and the - // possible list grows with the range - it reverted two healthy crystals and rescued none. - // - // 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 - // completeness and CC1/2 in different groups are not the same measurement. The search merge is - // P1 on both sides, over the full resolution range and before any correction surface is fitted - // - the one merge the two passes make in the same terms. + // 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) { - const double compl1 = pass1.search_merge_completeness; - const double compl2 = pass2.search_merge_completeness; - const bool compl_measured = pass1.search_merge_completeness_measured - && pass2.search_merge_completeness_measured; - // Say so rather than printing a zero that reads as a measurement. - const auto compl_text = [&] { - return compl_measured ? fmt::format("completeness {:.1f}% vs {:.1f}%", compl2, compl1) - // "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("completeness not compared"); - }; - const double cc1 = pass1.search_merge_cc_half; - const double cc2 = pass2.search_merge_cc_half; - constexpr double MAX_CREDIBLE_COMPLETENESS = 100.5; // rounding headroom, nothing more - - // Third arm: the axial rows. CC1/2 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 CC1/2 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, which is exactly what a wedge-losing - // pass looks like from CC1/2 - better. 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 only decides where CC1/2 does not: a pass that wins on CC1/2 by more than the loss - // the arm above already tolerates has said something about the bulk of the data that a few - // dozen reflections cannot answer, and is kept. And it needs pass 1 to have measured enough - // of the rows for the comparison to mean anything - below eight there is no zone a screw - // could be claimed on either side, so the difference decides nothing. - constexpr double MIN_AXIAL_ROWS_KEPT = 0.75; - constexpr int64_t MIN_AXIAL_ROWS_TO_COMPARE = 8; - const int64_t axial1 = pass1.search_merge_axial_reflections; - const int64_t axial2 = pass2.search_merge_axial_reflections; - const bool lost_axial_rows = axial1 >= MIN_AXIAL_ROWS_TO_COMPARE - && axial2 < MIN_AXIAL_ROWS_KEPT * static_cast(axial1) - && cc2 < cc1 + MAX_CC_HALF_LOSS; - - if (pass2.lattice_conflicts_with_fixed_sg - || (compl_measured && compl2 > MAX_CREDIBLE_COMPLETENESS) - || cc2 < cc1 - MAX_CC_HALF_LOSS || lost_axial_rows) { - const std::string worse = (cc2 < cc1 - MAX_CC_HALF_LOSS || !lost_axial_rows) - ? fmt::format("{}, CC1/2 before corrections {:.3f} vs {:.3f}", - compl_text(), cc2, cc1) - : fmt::format("it holds {} of the {} low-order axial reflections the " - "header-geometry pass measured, and its CC1/2 before corrections " - "{:.3f} vs {:.3f} does not say it is better", - axial2, axial1, cc2, cc1); + 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); @@ -1903,6 +1925,9 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { experiment_.Goniometer(restored); } 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(); auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); redo.post_refine = pass1.post_refine; redo.pass_number = pass2.pass_count + 1; @@ -1914,7 +1939,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { if (pass2.pass_decision.empty()) pass2.pass_decision = fmt::format( "post-refined geometry adopted ({}, CC1/2 before corrections {:.3f} vs {:.3f})", - compl_text(), cc2, cc1); + SearchMergeCompletenessText(pass1, pass2), pass2.search_merge_cc_half, + pass1.search_merge_cc_half); } if (pass2.pass_decision.empty()) pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_) @@ -2105,6 +2131,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const bool calibration_spots = calibration && config_.calibration_method == CalibrationMethod::Spots; const bool per_image_analysis = full || calibration; const bool write_files = write_output && !config_.output_prefix.empty(); + // Set below the first merge, once the pass can tell that RunAllPasses is going to throw it away and + // run again - see there for what that turns off and why it is knowable at that point. + bool superseded = false; // Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization, // space group, unit cell, ...) are configured on experiment_ by the caller. @@ -5161,12 +5190,36 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } } + // A pass RunAllPasses is about to throw away and run again. Its report is never printed and + // its files are overwritten seconds later, so everything below that only fills those in would + // be spent on an answer nobody reads. Every trigger of a re-run is settled by here - the + // integrator's background-starvation count (measured at integration), the fixed-group lattice + // conflict (at indexing) and the search merge the quality guard reads (at the first merge) - so + // a pass can tell that it is being superseded, which RunAllPasses cannot until it returns. + // Only the refined passes ask: quality_guard_pass1_ holds the header-geometry pass they are + // judged against and is set only while they run. The supercell guard is deliberately not asked + // here - it compares this pass's cell against pass 1's rather than reading this pass's own + // numbers - so a pass it re-runs still writes a report first. + // + // Nothing below this point is read by a later pass: the analyses are report-only (see each), + // the files are rewritten by the pass that supersedes this one, and the two places that touch + // experiment_ - the model frame and the P1 merge - are undone before the pass ends, the first + // by the SetSpaceGroup at the top of RunPipeline and the second on the spot. + superseded = + !cancelled_ && quality_guard_pass1_ && quality_guard_pass1_->has_merge_statistics + && ((bragg_adaptive_ && bkg_starved_fraction_ + && *bkg_starved_fraction_ > spot_width::BKG_STARVED_MAX_FRACTION) + || !RefinedPassIsWorse(*quality_guard_pass1_, result).empty()); + if (superseded) + logger.Info("This pass is going to be re-run, so it makes no report and writes no files."); + const auto &twin_sg_opt = experiment_.GetGemmiSpaceGroup(); const gemmi::SpaceGroup *twin_sg = twin_sg_opt ? &*twin_sg_opt : nullptr; - // Not on the geometry pre-pass: the analysis goes into that pass's statistics text and its - // written reflections, and neither survives the run. The promotion flag below is a different - // thing - it is what the SEARCH did, the second pass reads it, and it is set either way. - if (!geometry_prepass) { + // Not on the geometry pre-pass, nor on a superseded one: the analysis goes into that pass's + // statistics text and its written reflections, and neither survives the run. The promotion flag + // below is a different thing - it is what the SEARCH did, the second pass reads it, and it is + // set either way. + if (!geometry_prepass && !superseded) { // Translational pseudo-symmetry first: its verdict is an input to the twinning analysis, // which must not read its own L-test when a non-half-integer pseudo-translation has biased // it. Reporting only - no reflection and no merge depends on the outcome. @@ -5306,8 +5359,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log, none - // of which the geometry pre-pass produces. - if (!geometry_prepass) { + // of which the geometry pre-pass or a superseded pass produces. + if (!geometry_prepass && !superseded) { const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged); sm.statistics.wilson_b = wilson.b; sm.statistics.wilson_b_correlation = wilson.correlation; @@ -5385,8 +5438,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // sweep (the per-image scaling the rotation merge already fits, binned by frame = dose). A // per-image CC that falls, and/or a mosaicity that rises, with frame number is the classic // radiation-damage signature - a data-quality-vs-dose read complementary to the fitted decay - // correction. The full per-image table is written to _scaling.txt for detail. - if (experiment_.IsRotationIndexing()) { + // correction. The full per-image table is written to _scaling.txt for detail. Both this + // and the sweep-quality report below are text, so a superseded pass makes neither. + if (experiment_.IsRotationIndexing() && !superseded) { const auto &outs = indexer->GetIntegrationOutcome(); const int nf = static_cast(outs.size()); constexpr int nb = 10; @@ -5477,7 +5531,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // with no reference MTZ - the alternative indexing. Both are relabelings of the same // measurements, so the statistics computed above are unaffected; what changes is that the file, // the R-factors and the maps then all describe one indexing instead of two. - if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass) { + if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass + && !superseded) { phase("Validating against model"); const auto data_sg = experiment_.GetGemmiSpaceGroup(); // With a reference MTZ the merohedral indexing was already resolved against it (rotation @@ -5535,12 +5590,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b experiment_.GetSpaceGroupOrP1(), config_.output_prefix, logger); } - // Not on the geometry pre-pass. Pass 1 exists to choose the space group and post-refine the - // geometry; pass 2 remakes these files seconds later at the refined geometry, and that is the - // answer anyone reads. Writing them twice costs 0.6 s of a 15 s run - a fifth of the merged - // mmCIF on a large crystal - for a file that is superseded before the run ends. The pass-2 - // quality guard is unaffected: has_merge_statistics is set well above this, at the merge. - if (result.consensus_cell && write_files && config_.write_merged && !geometry_prepass) { + // Not on the geometry pre-pass, and not on a pass that is going to be re-run. Pass 1 exists to + // choose the space group and post-refine the geometry; the pass after it remakes these files + // seconds later at the refined geometry, and that is the answer anyone reads. Writing them + // twice costs 0.6 s of a 15 s run - a fifth of the merged mmCIF on a large crystal - for a file + // that is superseded before the run ends. The quality guard is unaffected: it reads the search + // merge, which is made well above this. + if (result.consensus_cell && write_files && config_.write_merged && !geometry_prepass + && !superseded) { phase("Writing reflections"); const ErrorModelReport em_report{ result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", @@ -5683,7 +5740,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // ones: the partiality and the per-image scale are left for the reading program, which fits a // scale model of its own. write_files, not write_output - without an output prefix there is // nowhere to put it, and the name would come out as a bare "_unmerged.mtz". - if (full && !cancelled_ && write_files && !geometry_prepass && result.consensus_cell) { + if (full && !cancelled_ && write_files && !geometry_prepass && !superseded && result.consensus_cell) { if (config_.export_unmerged) { if (observer) observer->OnPhase("Writing unmerged reflections"); const std::string path = config_.output_prefix + "_unmerged.mtz"; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index cdb07952e..1da6a0787 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -476,6 +476,11 @@ class Rugnux { // reflections left their background ring below six clean pixels. Measured, not predicted - it is // the only thing that tells a radius the pattern can take from one it cannot (see RunAllPasses). std::optional bkg_starved_fraction_; + // The header-geometry pass, for as long as the refined passes that are judged against it run (see + // the quality guard in RunAllPasses). Held here so a refined pass can ask the guard's own question + // of itself as soon as its search merge is in, which is where it learns that its report and its + // files are about to be superseded. + std::optional quality_guard_pass1_; // Pre-scan: read a spread sample of frames and take three things off them - the shadow of the beam // stop and its holder, added to the pixel mask (config_.detect_beam_stop), the beam centre