diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index f4550b8e4..842906684 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -6281,6 +6281,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } }; + // The geometry pre-pass integrates only the first prepass_fraction of the sweep (pass 2 reads it all). + const int prepass_end = start_image + + static_cast(std::lround((end_image - start_image) * config_.prepass_fraction)); + const int prepass_images = (prepass_end - start_image + config_.stride - 1) / config_.stride; + const int loop_end = geometry_prepass ? prepass_end : end_image; + auto full_worker = [&]() { pin_gpu(); // round-robin per worker thread; must precede engine construction // Before the analysis engine, which page-locks these bytes for its uploads and unregisters @@ -6291,9 +6297,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b AzimuthalIntegrationProfile profile(mapping); HarmonicEvidence worker_harmonic; std::vector worker_offsets; - const int loop_end = geometry_prepass - ? start_image + static_cast(std::lround((end_image - start_image) * config_.prepass_fraction)) - : end_image; while (!cancelled_) { const int ordinal = next_ordinal.fetch_add(1); @@ -6424,6 +6427,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b Note(fmt::format("{} images, {:.0f} per second", result.images_processed, result.images_processed / std::max(1e-6, result.image_loop_time_s))); + // A pre-pass that stopped at loop_end is a sweep of the frames it read, and that is all the merge and + // the post-refinement are handed. Left in, the empty tail is a stretch of frames every per-frame + // smoothing in the merge extrapolates over: measured, the fulls' log-scale curve, carried straight + // across 900 empty frames, fell further every round until the scale there was zero, and the + // cross-validation of its smoothness then took every scale in the sweep to NaN. + if (geometry_prepass) + indexer->GetIntegrationOutcome().resize(prepass_images); + // The footprint of the reflection rather than of the spot: the pre-scan's widths plus where the // spots sit against their predictions (SpotFootprint.h). Like the widths it acts only where the // reflection outgrows the r1 disk, and like them it is held back for the canonical pass. @@ -6984,14 +6995,27 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto &rot_ss = experiment_.GetScalingSettings(); const bool is_rotation = experiment_.IsRotationIndexing(); // rotation indexing -> rotation scaling/merge std::optional rsm; + std::optional guard_merge; // what the quality guard reads, where not the search merge std::optional prepass_postrefine_obs; // The rotation geometry post-refinement (see its call sites below for what it is for). const auto post_refine_geometry = [&] { const auto &outcomes = indexer->GetIntegrationOutcome(); if (geometry_prepass) { prepass_mosaicity_.assign(outcomes.size(), NAN); + std::vector fitted; for (size_t o = 0; o < outcomes.size(); ++o) - if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg; + if (outcomes[o].mosaicity_deg) { + prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg; + fitted.push_back(*outcomes[o].mosaicity_deg); + } + // A pre-pass of part of the sweep has fitted nothing for the frames after it, and they are + // predicted at the median of what it did fit, as the merge fills its own frames without a + // fit. Left at the frame's own estimate alone, the second half of a small-molecule sweep + // predicted a fifth fewer reflections and merged at R_meas 7.2 % against 5.2 %. + if (!fitted.empty() && outcomes.size() < static_cast(images_to_process)) { + std::nth_element(fitted.begin(), fitted.begin() + fitted.size() / 2, fitted.end()); + prepass_mosaicity_.resize(images_to_process, fitted[fitted.size() / 2]); + } } if (result.consensus_cell.has_value()) { if (const auto rot = indexer->FinalizeRotationIndexing(); rot && rot->axis) { @@ -7118,6 +7142,38 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Rotation scaling/merging (RotationScaleMerge) does not support " "wedge refinement"); + // The quality guard judges a refined pass against the pre-pass on the search merge, and a + // pre-pass that integrated only the first part of the sweep merged only that part. Against + // it a merge of the whole sweep holds more reflections, more signal and more of the axial + // rows for no reason but the frames, and the guard could never find it worse. So the guard + // reads this pass's own integration of the same frames, merged the way the pre-pass merged + // them. The merge writes its per-frame scale back, so it gets copies of the outcomes and + // only borrows the reflections. + if (quality_guard_pass1_ && prepass_images < images_to_process) { + logger.Info("Two-pass: the quality guard reads this pass's merge of the first {} images, " + "the ones the pre-pass integrated", prepass_images); + auto &outcomes = indexer->GetIntegrationOutcome(); + std::vector first_part(prepass_images); + for (int i = 0; i < prepass_images; ++i) { + ReflectionVector reflections = std::move(outcomes[i].reflections); + first_part[i] = outcomes[i]; + first_part[i].reflections = std::move(reflections); + } + { + RotationScaleMerge first_part_rsm(experiment_, first_part, result.consensus_cell, + static_cast(config_.scaling_iter), config_.nthreads, + logger, ""); + first_part_rsm.Ingest(); + first_part_rsm.SetExportScaledFulls(false); + first_part_rsm.SetFrenchWilson(false); + auto r = first_part_rsm.Run(!experiment_.GetGemmiSpaceGroup().has_value(), + /*full_stats=*/false, /*measure_cc_before_corrections=*/true); + guard_merge = ScaleMergeResult{std::move(r.merged), std::move(r.statistics), + r.cc_half_before_corrections}; + } + for (int i = 0; i < prepass_images; ++i) + outcomes[i].reflections = std::move(first_part[i].reflections); + } // A reference MTZ is allowed for rotation: it fixes the space group / cell (on the CLI) and // resolves the indexing ambiguity (below), but is NOT used to scale - the rotation merge stays // self-consistent. @@ -7279,15 +7335,17 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // The two-pass quality guard judges the second pass against the first on THIS merge, which both // passes make in the same terms - P1 (or the group the user fixed), full range, no correction - // surfaces - where their final merges need not even be in the same group. See RunAllPasses. + // surfaces - where their final merges need not even be in the same group. See RunAllPasses. Where + // the pre-pass integrated only part of the sweep, it is the same merge of those frames (above). { - const auto &o = sm.statistics.overall; + const ScaleMergeResult &counted = guard_merge ? *guard_merge : sm; + const auto &o = counted.statistics.overall; result.search_merge_completeness_measured = o.possible_unique_reflections > 0; result.search_merge_completeness = result.search_merge_completeness_measured ? 100.0 * o.unique_reflections / o.possible_unique_reflections : 0.0; - result.search_merge_cc_half = sm.cc_half_before_corrections; - if (!sm.statistics.shells.empty()) - result.search_merge_inner_cc_half = sm.statistics.shells.front().cc_half; + result.search_merge_cc_half = counted.cc_half_before_corrections; + if (!counted.statistics.shells.empty()) + result.search_merge_inner_cc_half = counted.statistics.shells.front().cc_half; // ... and the reflections it holds on the principal axial rows, which the guard needs // separately because CC1/2 cannot express them (see RunAllPasses). Only the low-order // part is counted: a screw is read off the strong start of its row - the evidence gate @@ -7316,7 +7374,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b }; constexpr float AXIAL_ROW_LOW_ORDER_D_A = 4.0f; result.search_merge_axial_reflections = std::count_if( - sm.merged.begin(), sm.merged.end(), [&](const MergedReflection &r) { + counted.merged.begin(), counted.merged.end(), [&](const MergedReflection &r) { return ((r.h == 0) + (r.k == 0) + (r.l == 0)) == 2 && r.d > AXIAL_ROW_LOW_ORDER_D_A && centring_allows(r); }); @@ -7324,14 +7382,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // what the guard decides on. constexpr float STRONG_I_OVER_SIGMA = 2.0f; result.search_merge_strong_reflections = std::count_if( - sm.merged.begin(), sm.merged.end(), [&](const MergedReflection &r) { + counted.merged.begin(), counted.merged.end(), [&](const MergedReflection &r) { return r.sigma > 0 && r.I / r.sigma >= STRONG_I_OVER_SIGMA && centring_allows(r); }); result.search_merge_negative_reflections = std::count_if( - sm.merged.begin(), sm.merged.end(), [&](const MergedReflection &r) { + counted.merged.begin(), counted.merged.end(), [&](const MergedReflection &r) { return r.sigma > 0 && r.I / r.sigma <= -STRONG_I_OVER_SIGMA && centring_allows(r); }); - result.search_merge_reflections = std::count_if(sm.merged.begin(), sm.merged.end(), + result.search_merge_reflections = std::count_if(counted.merged.begin(), counted.merged.end(), centring_allows); }