diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index c4cdca5c..0134e2bb 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,9 @@ # Changelog ## 1.0.0 ### Unreleased +* **rugnux: the first pass of the rotation two-pass no longer merges, reports or writes anything.** It exists to measure the detector geometry, the goniometer rotation scale, the mosaicity and the space group, and the second pass makes the merged result again at the refined geometry - so the first pass now stops once it has those. **The `_01_*` files are no longer written.** A battery of 24 rotation crystals drops from 6m52s to 6m17s and no merged result changes. +* rugnux: every run ends with a `Time by phase` table - wall time and mean cores busy for each phase of the whole run, both passes together - so a slow dataset can be diagnosed from its own log without a profiler. +* rugnux: **fixed** the `Per-image cost` line, which divided the per-stage means by the thread count while the image loop runs a few workers per GPU. On a 2-GPU machine at `-N 48` every stage was reported 6x too small. It now divides by, and names, the loop's own worker count. * **rugnux: scaling and merging are much faster on crystals that integrate far beyond the resolution they merge at.** Observations outside the scaling resolution range are dropped as they are ingested instead of being scaled, combined and post-refined first, and the geometry post-refinement fits on a bounded sample of them. On a large-cell dataset merging at less than half the resolution its detector reaches, a run drops from 68 s to 37 s; merged statistics are unchanged. * rugnux: the detector-frame modulation correction is fitted on a grid spanning the detector rather than the reflections that happen to be present, so whether it is applied no longer depends on how far integration reached. On a crystal where it was being refused, merged R_meas improves by 4 percentage points. * rugnux: the first-pass rotation indexing finds its spots on every worker rather than one, which is worth most on large detectors - on a 16M-pixel dataset that phase drops by 44% and the whole run by 9%. The lattice it picks is unchanged. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 0605fd30..5f725eec 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1,6 +1,10 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#ifdef __linux__ +#include +#endif + #include "Rugnux.h" #include "ModelValidation.h" @@ -331,8 +335,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // Say which of the two this is reading for. It runs before the first image of the run proper is // processed, so without this the job sits at its initial status for the whole pre-scan - which is // tens of seconds once the beam centre asks for its own frames on top of the projection. - if (observer) - observer->OnPhase(want_shadow ? (want_beam_center ? "Beam stop and beam centre" : "Beam stop") + MarkPhase(observer, want_shadow ? (want_beam_center ? "Beam stop and beam centre" : "Beam stop") : "Beam centre"); ShadowFinder finder(experiment_, pixel_mask_); @@ -518,8 +521,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru frame_angle_deg.size(), beam_center_spots.size(), extra.size()); // Its own phase: this reads several times the frames the pass above did, and is the one // step of the pre-scan a user is likely to sit through wondering what is happening. - if (observer) - observer->OnPhase("Beam centre (reading more frames)"); + MarkPhase(observer, "Beam centre (reading more frames)"); // Read on workers as above - this pass is several times the size of the first one - and // append in `extra` order so the pool does not depend on how they interleaved. std::vector> extra_spots(extra.size()); @@ -606,8 +608,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ } const int refine_frames = std::max(1, config_.refine_geometry.value()); - if (observer) - observer->OnPhase("Geometry refinement (first pass)"); + MarkPhase(observer, "Geometry refinement (first pass)"); const auto dataset = reader_.GetDataset(); @@ -775,7 +776,6 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { if (config_.mode == ProcessMode::FullAnalysis && config_.rotation_postrefine_geometry && experiment_.IsRotationIndexing()) { Logger logger("Rugnux"); - const std::string base_prefix = config_.output_prefix; const auto gonio_snapshot = experiment_.GetGoniometer(); prepass_detector_geometry_.reset(); prepass_rotation_scale_.reset(); @@ -784,15 +784,12 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { prepass_merge_sg_.reset(); prepass_promoted_point_group_ = false; - 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); + logger.Info("Rotation two-pass geometry post-refinement: first pass (header geometry), measuring " + "the geometry the second pass will integrate at"); + auto pass1 = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/true); pass1.pass_number = 1; pass1.pass_count = 2; - if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } + if (cancelled_) return pass1; if (gonio_snapshot) experiment_.Goniometer(*gonio_snapshot); // undo the pre-pass goniometer shift // Apply the post-refined detector geometry for the second pass, keeping the header geometry so @@ -834,11 +831,10 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // the centring and pass 2 dropped it, and the run merged in P1. experiment_.SpaceGroupNumber(std::nullopt); prepass_merge_sg_ = pass1.space_group_number; - prepass_promoted_point_group_ = pass1.twinning.laue_class_was_chosen_by_promotion; + prepass_promoted_point_group_ = pass1.promoted_point_group; 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) + config_.output_prefix); auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); // 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. @@ -919,7 +915,6 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { if (prepass_rotation_scale_ && gonio_snapshot) experiment_.Goniometer(*gonio_snapshot); // and the header rotation angles with it experiment_.SpaceGroupNumber(std::nullopt); - config_.output_prefix = base_prefix; auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); if (!redo.space_group_search.has_value()) redo.space_group_search = pass1.space_group_search; @@ -1000,6 +995,50 @@ namespace { } } +namespace { + // CPU time this process has used across all its threads, in seconds. Compared against wall time it + // says how many cores a phase kept busy, which is the difference between a phase that is slow + // because it has work to do and one that is slow because it is running on a single thread. + double process_cpu_seconds() { +#ifdef __linux__ + rusage ru{}; + if (getrusage(RUSAGE_SELF, &ru) != 0) + return 0.0; + const auto secs = [](const timeval &t) { return t.tv_sec + t.tv_usec * 1e-6; }; + return secs(ru.ru_utime) + secs(ru.ru_stime); +#else + return 0.0; +#endif + } +} + +void Rugnux::ClosePhase() { + if (phase_current_.empty()) + return; + const double wall = std::chrono::duration(std::chrono::steady_clock::now() - phase_start_).count(); + const double cpu = process_cpu_seconds() - phase_start_cpu_s_; + // A phase entered again - the two passes of a rotation run enter most of them twice - adds to the + // entry it already has, so the report has one line per phase rather than one per visit. + auto it = std::find_if(phase_timings_.begin(), phase_timings_.end(), + [this](const PhaseTiming &p) { return p.name == phase_current_; }); + if (it == phase_timings_.end()) + phase_timings_.push_back({phase_current_, wall, cpu}); + else { + it->wall_s += wall; + it->cpu_s += cpu; + } + phase_current_.clear(); +} + +void Rugnux::MarkPhase(RugnuxObserver *observer, const std::string &phase) { + ClosePhase(); + phase_current_ = phase; + phase_start_ = std::chrono::steady_clock::now(); + phase_start_cpu_s_ = process_cpu_seconds(); + if (observer) + observer->OnPhase(phase); +} + ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) { Logger logger("Rugnux"); ProcessResult result; @@ -1123,8 +1162,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b : calibration ? "powder calibration" : "azimuthal integration"; logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]", images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name); - if (observer) - observer->OnPhase(full ? "Full analysis" : calibration ? "Powder calibration" : "Azimuthal integration"); + MarkPhase(observer, full ? "Setting up" : calibration ? "Powder calibration" : "Azimuthal integration"); // Full-analysis shared engines. std::unique_ptr indexer_pool; @@ -1166,8 +1204,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b indexer->ForceRotationIndexerLattice(*config_.forced_rotation_lattice); logger.Info("Rotation indexer lattice forced externally - skipping first pass"); } else if (full && config_.rotation_indexing && config_.two_pass_rotation) { - if (observer) - observer->OnPhase("Rotation indexing (first pass)"); + MarkPhase(observer, "Rotation indexing (first pass)"); // Mid-exposure goniometer angle for an image ordinal (matches the per-image path). auto rot_angle = [&](int ordinal) -> std::optional { @@ -1755,8 +1792,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } }; - if (observer) - observer->OnPhase("Processing images"); + MarkPhase(observer, "Processing images"); std::function worker = per_image_analysis ? std::function(full_worker) : std::function(azint_worker); @@ -1769,6 +1805,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const int loop_workers = get_gpu_count() > 0 ? std::min(config_.nthreads, std::max(8, 4 * get_gpu_count())) : config_.nthreads; + // The per-stage means below are per WORKER, so anything reporting them has to divide by this and + // not by the thread count - on a machine with GPUs they are not the same number. + result.loop_workers = loop_workers; std::vector > futures; futures.reserve(loop_workers); const auto image_loop_start = std::chrono::steady_clock::now(); @@ -1803,8 +1842,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and // writes the PONI file; nothing here is stored in the _process.h5. if (calibration && !cancelled_) { - if (observer) - observer->OnPhase("Powder calibration"); + MarkPhase(observer, "Powder calibration"); if (calibration_spots) { logger.Info("Powder calibration from {} pooled spots", calibration_spot_list.size()); result.calibration = CalibrateFromSpots(calibration_spot_list, @@ -1865,9 +1903,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b && (config_.run_scaling || !config_.reference_data.empty())) { // Scaling/merging is a long post-pass; report each sub-step as a phase so the GUI progress // bar reflects what is happening instead of freezing on one label. - auto phase = [&](const std::string &p) { - if (observer) observer->OnPhase(p); - }; + auto phase = [&](const std::string &p) { MarkPhase(observer, p); }; phase("Scaling and merging"); // Ice-ring handling (--detect-ice-rings): flag reflections sitting on a hexagonal-ice powder @@ -2457,10 +2493,15 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b experiment_.SpaceGroupNumber(sg.number); end_msg.space_group_number = sg.number; result.space_group_number = sg.number; - phase("Re-scaling in space group " + sg.short_name()); - sm = scale_and_merge(sg.short_name(), false); + // The pre-pass wanted the GROUP, which is now determined. Merging in it is what the + // second pass does, at the refined geometry, and that merge is the one anybody reads. + if (!geometry_prepass) { + phase("Re-scaling in space group " + sg.short_name()); + sm = scale_and_merge(sg.short_name(), false); + } } result.space_group_search = sg_search; + result.promoted_point_group = promoted_point_group; } else { // A space group was fixed by the user; surface it so the viewer/CLI can still show it. result.space_group_number = experiment_.GetSpaceGroupNumber(); @@ -2470,287 +2511,294 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // second pass REUSES it instead of re-searching at the slightly changed geometry (which could flip a // borderline determination and mix the two passes). The lattice is deliberately NOT forced - the second // pass re-indexes at the refined geometry so the CELL comes out self-consistent with it (forcing the - // pre-pass lattice would pin the nominal cell and undo the geometry refinement). Unlike a throwaway - // pre-pass, this one now CONTINUES to the write below so its own (header-geometry) result is saved - // as the first pass. + // pre-pass lattice would pin the nominal cell and undo the geometry refinement). + // + // And that is everything the pre-pass exists to measure: the post-refined detector geometry, the + // goniometer rotation scale, the frame-order-smoothed mosaicity, the space group and (for the + // supercell guard) its indexing result. Everything below produces the merged intensities and the + // report and files that go with them, and the second pass makes all of it again at the refined + // geometry - so for the pre-pass it is a merge nobody reads. Measured over the rotation test set, + // running it anyway was 10% of the battery. - // Reference-based indexing-ambiguity resolution (rotation). When a reference MTZ is supplied and - // the crystal has an indexing ambiguity (merohedral: lattice symmetry higher than the Laue group), - // pick the reindexing that best correlates with the reference intensities and re-merge in it. The - // twin-law reindex is metric-preserving, so only the hkl labels change (the cell is unchanged). - // Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here. - if (rsm && !config_.reference_data.empty() && result.consensus_cell - && experiment_.GetSpaceGroupNumber().has_value()) { - const int sg_num = static_cast(*experiment_.GetSpaceGroupNumber()); - const auto choice = ChooseReindex( - sm.merged, *result.consensus_cell, sg_num, - [&](const std::vector &m) { - return ReferenceIntensityCC(m, config_.reference_data, sg_num); - }); - if (!choice.is_identity) { - logger.Info("Reference: resolved indexing ambiguity by reindexing to match the reference " - "(CC {:.3f}, vs {:.3f} unchanged)", choice.score, choice.identity_score); - for (auto &io : indexer->GetIntegrationOutcome()) - for (auto &r : io.reflections) { - const gemmi::Op::Miller h = choice.op.apply_to_hkl({{r.h, r.k, r.l}}); - r.h = h[0]; r.k = h[1]; r.l = h[2]; - } - rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell, - static_cast(config_.scaling_iter), - config_.nthreads, logger, config_.observation_dump_path); - rsm->Ingest(); - phase("Re-merging in the reference indexing"); - sm = scale_and_merge(experiment_.GetGemmiSpaceGroup()->short_name(), false); + if (!geometry_prepass) { + // Reference-based indexing-ambiguity resolution (rotation). When a reference MTZ is supplied and + // the crystal has an indexing ambiguity (merohedral: lattice symmetry higher than the Laue group), + // pick the reindexing that best correlates with the reference intensities and re-merge in it. The + // twin-law reindex is metric-preserving, so only the hkl labels change (the cell is unchanged). + // Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here. + if (rsm && !config_.reference_data.empty() && result.consensus_cell + && experiment_.GetSpaceGroupNumber().has_value()) { + const int sg_num = static_cast(*experiment_.GetSpaceGroupNumber()); + const auto choice = ChooseReindex( + sm.merged, *result.consensus_cell, sg_num, + [&](const std::vector &m) { + return ReferenceIntensityCC(m, config_.reference_data, sg_num); + }); + if (!choice.is_identity) { + logger.Info("Reference: resolved indexing ambiguity by reindexing to match the reference " + "(CC {:.3f}, vs {:.3f} unchanged)", choice.score, choice.identity_score); + for (auto &io : indexer->GetIntegrationOutcome()) + for (auto &r : io.reflections) { + const gemmi::Op::Miller h = choice.op.apply_to_hkl({{r.h, r.k, r.l}}); + r.h = h[0]; r.k = h[1]; r.l = h[2]; + } + rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell, + static_cast(config_.scaling_iter), + config_.nthreads, logger, config_.observation_dump_path); + rsm->Ingest(); + phase("Re-merging in the reference indexing"); + sm = scale_and_merge(experiment_.GetGemmiSpaceGroup()->short_name(), false); + } } - } - const auto twin_sg_number = experiment_.GetSpaceGroupNumber(); - const gemmi::SpaceGroup *twin_sg = twin_sg_number - ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; - result.twinning = AnalyzeTwinning(sm.merged, twin_sg); - // Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the - // search itself made, so the text cannot claim "no twin law exists" on its own say-so. - result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group; - stats_text << TwinningAnalysisToText(result.twinning) << "\n"; + const auto twin_sg_number = experiment_.GetSpaceGroupNumber(); + const gemmi::SpaceGroup *twin_sg = twin_sg_number + ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; + result.twinning = AnalyzeTwinning(sm.merged, twin_sg); + // Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the + // search itself made, so the text cannot claim "no twin law exists" on its own say-so. + result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group; + stats_text << TwinningAnalysisToText(result.twinning) << "\n"; - // Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the - // crystal was mounted, which the user can change), warned about when the two nearly coincide. - if (experiment_.IsRotationIndexing() && twin_sg && end_msg.rotation_lattice.has_value()) { - if (const auto gonio = experiment_.GetGoniometer()) { - const auto closest = ClosestSymmetryAxisToSpindle(*twin_sg, *end_msg.rotation_lattice, - gonio->GetAxis()); - if (closest.has_value()) { - // Practical bound, not a derived one: the blind cone's half-angle is the maximum - // Bragg angle (~15 deg for 2 A data at 1 A), and measured on this battery a 13.6 deg - // case shows the loss while a 16.2 deg one is 99.7% complete. - constexpr double WARN_DEG = 15.0; - const auto [angle, order] = *closest; - if (angle < WARN_DEG) { + // Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the + // crystal was mounted, which the user can change), warned about when the two nearly coincide. + if (experiment_.IsRotationIndexing() && twin_sg && end_msg.rotation_lattice.has_value()) { + if (const auto gonio = experiment_.GetGoniometer()) { + const auto closest = ClosestSymmetryAxisToSpindle(*twin_sg, *end_msg.rotation_lattice, + gonio->GetAxis()); + if (closest.has_value()) { + // Practical bound, not a derived one: the blind cone's half-angle is the maximum + // Bragg angle (~15 deg for 2 A data at 1 A), and measured on this battery a 13.6 deg + // case shows the loss while a 16.2 deg one is 99.7% complete. + constexpr double WARN_DEG = 15.0; + const auto [angle, order] = *closest; + if (angle < WARN_DEG) { + const std::string msg = fmt::format( + "The crystal's {}-fold axis is only {:.1f} deg from the spindle. A rotation sweep " + "never records the reflections whose reciprocal vector lies within the Bragg angle " + "of the spindle, and symmetry normally supplies them from an equivalent elsewhere; " + "it cannot here, because that axis maps the blind region onto itself. Those " + "reflections stay missing however long the sweep runs. The loss is confined to " + "that cone rather than spread over the data, so overall completeness may still " + "look reasonable - check it near the spindle direction. A second sweep on a " + "different axis, or re-mounting, recovers them.", + order, angle); + logger.Warning("{}", msg); + stats_text << " !! " << msg << "\n\n"; + result.warnings.push_back(msg); + } else { + stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at " + << std::fixed << std::setprecision(1) << angle << " deg\n\n"; + } + } + } + } + + // Indexing-ambiguity (alternative-indexing) advisory. When the lattice metric symmetry exceeds + // the Laue symmetry the crystal can be validly indexed in several hands related by twin-law + // operators. For an obvious merohedral case (P3/P4/P6...) users expect this; but a PSEUDO-merohedral + // metric (e.g. a C2 crystal whose beta makes it pseudo-F-orthorhombic) is easy to miss, so surface it. The + // reindex operator is applied directly to (h,k,l), so its triplet reads as an h,k,l transform. + if (result.consensus_cell && twin_sg_number) { + const auto twin_ops = ReindexAmbiguityOperators(*result.consensus_cell, + static_cast(*twin_sg_number), 2.0); + if (!twin_ops.empty()) { + std::string laws; + for (const auto &op : twin_ops) { + std::string t = op.triplet(); + std::replace(t.begin(), t.end(), 'x', 'h'); + std::replace(t.begin(), t.end(), 'y', 'k'); + std::replace(t.begin(), t.end(), 'z', 'l'); + laws += (laws.empty() ? "" : " ; ") + t; + } + if (config_.reference_data.empty()) { const std::string msg = fmt::format( - "The crystal's {}-fold axis is only {:.1f} deg from the spindle. A rotation sweep " - "never records the reflections whose reciprocal vector lies within the Bragg angle " - "of the spindle, and symmetry normally supplies them from an equivalent elsewhere; " - "it cannot here, because that axis maps the blind region onto itself. Those " - "reflections stay missing however long the sweep runs. The loss is confined to " - "that cone rather than spread over the data, so overall completeness may still " - "look reasonable - check it near the spindle direction. A second sweep on a " - "different axis, or re-mounting, recovers them.", - order, angle); + "Indexing ambiguity: this cell / space group admits alternative indexing " + "(reindex operator(s): {}). {} are indexed in one hand at random, and rugnux can only " + "break this against an external reference. WITHOUT one the merge mixes the hands and " + "CC1/2 is degraded - supply a reference (-z reference.mtz or --model) to resolve it.", + laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals"); logger.Warning("{}", msg); stats_text << " !! " << msg << "\n\n"; result.warnings.push_back(msg); } else { - stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at " - << std::fixed << std::setprecision(1) << angle << " deg\n\n"; + logger.Info("Indexing ambiguity present (reindex operator(s): {}); resolved against the " + "supplied reference.", laws); } } } - } - // Indexing-ambiguity (alternative-indexing) advisory. When the lattice metric symmetry exceeds - // the Laue symmetry the crystal can be validly indexed in several hands related by twin-law - // operators. For an obvious merohedral case (P3/P4/P6...) users expect this; but a PSEUDO-merohedral - // metric (e.g. a C2 crystal whose beta makes it pseudo-F-orthorhombic) is easy to miss, so surface it. The - // reindex operator is applied directly to (h,k,l), so its triplet reads as an h,k,l transform. - if (result.consensus_cell && twin_sg_number) { - const auto twin_ops = ReindexAmbiguityOperators(*result.consensus_cell, - static_cast(*twin_sg_number), 2.0); - if (!twin_ops.empty()) { - std::string laws; - for (const auto &op : twin_ops) { - std::string t = op.triplet(); - std::replace(t.begin(), t.end(), 'x', 'h'); - std::replace(t.begin(), t.end(), 'y', 'k'); - std::replace(t.begin(), t.end(), 'z', 'l'); - laws += (laws.empty() ? "" : " ; ") + t; - } - if (config_.reference_data.empty()) { - const std::string msg = fmt::format( - "Indexing ambiguity: this cell / space group admits alternative indexing " - "(reindex operator(s): {}). {} are indexed in one hand at random, and rugnux can only " - "break this against an external reference. WITHOUT one the merge mixes the hands and " - "CC1/2 is degraded - supply a reference (-z reference.mtz or --model) to resolve it.", - laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals"); - logger.Warning("{}", msg); - stats_text << " !! " << msg << "\n\n"; - result.warnings.push_back(msg); - } else { - logger.Info("Indexing ambiguity present (reindex operator(s): {}); resolved against the " - "supplied reference.", laws); - } - } - } - - // 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. - { - const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged); - sm.statistics.wilson_b = wilson.b; - sm.statistics.wilson_b_correlation = wilson.correlation; - if (std::isfinite(wilson.b) && wilson.b > 0.0) - logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)", - wilson.b, wilson.correlation, wilson.n_shells); - } - - stats_text << sm.statistics; - result.merge_statistics_text = stats_text.str(); - result.has_merge_statistics = true; - result.merge_statistics = sm.statistics; - - // Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where - // the image is in no flagged range, plus the vocabulary the codes index. Only filled when the - // diagnostic ran, so absent datasets mean "not looked for" rather than "all clean". - if (sm.statistics.sweep_quality.measured) { - end_msg.sweep_quality.assign(end_msg.max_image_number, 0); - for (const auto &r : sm.statistics.sweep_quality.ranges) - for (int64_t i = std::max(0, r.first_image); - i <= r.last_image && i < static_cast(end_msg.sweep_quality.size()); ++i) - end_msg.sweep_quality[i] = static_cast(r.reason) + 1; - for (int r = 0; r <= static_cast(SweepQualityReason::RadiationDamage); ++r) - end_msg.sweep_quality_reasons.emplace_back( - SweepQualityReasonCode(static_cast(r))); - } - { - // Stride rather than take the head: the merged list is ordered by hkl, so the first N - // reflections are one corner of reciprocal space and would not show the intensity range. - constexpr size_t MAX_POINTS = 4000; - const size_t n = sm.merged.size(); - const size_t step = std::max(1, n / MAX_POINTS); - result.merged_i_sigma.reserve(std::min(n, MAX_POINTS + 1)); - for (size_t i = 0; i < n; i += step) - if (std::isfinite(sm.merged[i].I) && std::isfinite(sm.merged[i].sigma) - && sm.merged[i].sigma > 0.0f) - result.merged_i_sigma.emplace_back(sm.merged[i].I, sm.merged[i].sigma); - } - result.has_reference = !config_.reference_data.empty(); - - // Inherit the campaign's shared R-free test set from the reference MTZ (overriding the - // per-hkl hash that the merge assigned), so every dataset flags the same free reflections. - if (config_.reference_has_free_flags && !config_.reference_data.empty() && !sm.merged.empty()) { - const auto sg = experiment_.GetSpaceGroupNumber().value_or(1); - const size_t matched = ApplyReferenceFreeFlags(sm.merged, static_cast(sg), - config_.reference_data); - logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)", - matched, sm.merged.size()); - } - - // Radiation-damage report (rotation): the per-image scale CC-to-merge and mosaicity across the - // 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()) { - const auto &outs = indexer->GetIntegrationOutcome(); - const int nf = static_cast(outs.size()); - constexpr int nb = 10; - std::array cc_sum{}, mos_sum{}; - std::array cc_n{}, mos_n{}; - for (int f = 0; f < nf; ++f) { - const int b = nf > 1 ? std::min(nb - 1, f * nb / nf) : 0; - if (outs[f].image_scale_cc && std::isfinite(*outs[f].image_scale_cc)) { - cc_sum[b] += *outs[f].image_scale_cc; cc_n[b]++; - } - if (outs[f].mosaicity_deg && std::isfinite(*outs[f].mosaicity_deg)) { - mos_sum[b] += *outs[f].mosaicity_deg; mos_n[b]++; - } - } - int tot = 0; - for (int n : cc_n) tot += n; - if (tot >= 20) { - std::ostringstream os; - os << fmt::format("Radiation-damage report (per-image scale over {} frames):\n", nf); - os << " dose \n"; - for (int b = 0; b < nb; ++b) { - if (cc_n[b] == 0 && mos_n[b] == 0) continue; - os << fmt::format(" {:3d}-{:3d}% {:>7} {:>7}\n", b * 10, (b + 1) * 10, - cc_n[b] ? fmt::format("{:.3f}", cc_sum[b] / cc_n[b]) : std::string("-"), - mos_n[b] ? fmt::format("{:.3f}", mos_sum[b] / mos_n[b]) : std::string("-")); - } - if (cc_n[0] > 0 && cc_n[nb - 1] > 0) { - const double cc0 = cc_sum[0] / cc_n[0], cc9 = cc_sum[nb - 1] / cc_n[nb - 1]; - os << fmt::format(" => per-image CC to merge {:.3f} (first 10%) -> {:.3f} (last 10%){}", - cc0, cc9, - cc9 < cc0 - 0.05 ? " (falling: radiation damage / crystal decay)" : ""); - } - // Relative B-factor across the sweep (measured before any decay correction, against the - // low-dose start): the resolution-dependent complement to the CC/mosaicity read above. - // Radiation damage fades the high-resolution intensity, so only a POSITIVE change is dose. - // A curve no straight line describes gets no first->last number at all, and a batch whose - // data could not measure one prints "-" rather than a value. - if (!sm.statistics.radiation_damage_b_batch.empty()) { - if (std::isfinite(sm.statistics.radiation_damage_delta_b)) - os << fmt::format("\n => relative B-factor change over run = {:+.2f} A^2 (first->last){}", - sm.statistics.radiation_damage_delta_b, - sm.statistics.radiation_damage_delta_b > 5.0 - ? " (significant dose-dependent scaling: radiation damage)" : ""); - else - os << "\n => relative B-factor: the per-batch curve below is not a trend, so no " - "first->last number describes it - dose does not come back, so whatever moved " - "here was not dose; see the sweep-quality report"; - os << fmt::format("\n per-batch relative-B (A^2, {:.0f} deg/batch):", - sm.statistics.radiation_damage_batch_deg); - for (float bb : sm.statistics.radiation_damage_b_batch) - os << (std::isfinite(bb) ? fmt::format(" {:.1f}", bb) : std::string(" -")); - } - logger.Info("{}", os.str()); - result.radiation_damage_text = os.str(); + // 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. + { + const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged); + sm.statistics.wilson_b = wilson.b; + sm.statistics.wilson_b_correlation = wilson.correlation; + if (std::isfinite(wilson.b) && wilson.b > 0.0) + logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)", + wilson.b, wilson.correlation, wilson.n_shells); } - // Sweep-quality report: the stretches of the sweep over which the crystal delivered much - // less than the rest of the run, and what each one looks like. Nothing is excluded because - // of it - the frames still carry signal, and this is a message for the beamline, not a - // filter. Frame numbers are processed-image ordinals, inclusive, as in _image.dat. - const auto &sq = sm.statistics.sweep_quality; - if (sq.measured) { - std::ostringstream os; - os << fmt::format("Sweep quality over {:.0f} deg (per-image scale with the incident flux, " - "which varied {:.2f}x, already divided out):\n", - sq.sweep_deg, sq.flux_peak_to_trough); - if (sq.ranges.empty()) { - os << " no stretch of the sweep is materially worse than the run"; - } else { - os << " frames rotation diagnosis severity scale CC indexed\n"; - for (const auto &r : sq.ranges) - os << fmt::format(" {:<17s} {:6.1f} deg {:<16s} {:5.2f} {:5.2f} {:5.2f} {:4.0f}%\n", - fmt::format("{}-{}", r.first_image, r.last_image), r.rotation_deg, - SweepQualityReasonText(r.reason), r.severity, r.mean_relative_scale, - r.mean_relative_cc, 100.0 * r.indexed_fraction); - os << " => severity is the fraction of the run's typical diffracting power missing " - "over the range"; - if (sq.modulation_peak_to_trough >= 1.05f) - os << fmt::format("\n => once-per-revolution modulation of the per-image scale: " - "{:.1f}x peak to trough", sq.modulation_peak_to_trough); - } - logger.Info("{}", os.str()); + stats_text << sm.statistics; + result.merge_statistics_text = stats_text.str(); + result.has_merge_statistics = true; + result.merge_statistics = sm.statistics; + + // Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where + // the image is in no flagged range, plus the vocabulary the codes index. Only filled when the + // diagnostic ran, so absent datasets mean "not looked for" rather than "all clean". + if (sm.statistics.sweep_quality.measured) { + end_msg.sweep_quality.assign(end_msg.max_image_number, 0); + for (const auto &r : sm.statistics.sweep_quality.ranges) + for (int64_t i = std::max(0, r.first_image); + i <= r.last_image && i < static_cast(end_msg.sweep_quality.size()); ++i) + end_msg.sweep_quality[i] = static_cast(r.reason) + 1; + for (int r = 0; r <= static_cast(SweepQualityReason::RadiationDamage); ++r) + end_msg.sweep_quality_reasons.emplace_back( + SweepQualityReasonCode(static_cast(r))); } - } + { + // Stride rather than take the head: the merged list is ordered by hkl, so the first N + // reflections are one corner of reciprocal space and would not show the intensity range. + constexpr size_t MAX_POINTS = 4000; + const size_t n = sm.merged.size(); + const size_t step = std::max(1, n / MAX_POINTS); + result.merged_i_sigma.reserve(std::min(n, MAX_POINTS + 1)); + for (size_t i = 0; i < n; i += step) + if (std::isfinite(sm.merged[i].I) && std::isfinite(sm.merged[i].sigma) + && sm.merged[i].sigma > 0.0f) + result.merged_i_sigma.emplace_back(sm.merged[i].I, sm.merged[i].sigma); + } + result.has_reference = !config_.reference_data.empty(); - if (result.consensus_cell && write_files && config_.write_merged) { - phase("Writing reflections"); - const ErrorModelReport em_report{ - result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", - result.error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", result.error_model_isa_asymptotic) - : std::string(), - result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(), - result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()}; - WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, - em_report, result.twinning, config_.output_prefix); - // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS - // comparison. The offline self-scaling result is otherwise not exposed (process.h5's - // per-image arrays are only filled on the online per-image path). Sourced from the - // partials, which carry the first-pass per-image scale. - ScalingResult(indexer->GetIntegrationOutcome()).SaveToFile(config_.output_prefix); - } + // Inherit the campaign's shared R-free test set from the reference MTZ (overriding the + // per-hkl hash that the merge assigned), so every dataset flags the same free reflections. + if (config_.reference_has_free_flags && !config_.reference_data.empty() && !sm.merged.empty()) { + const auto sg = experiment_.GetSpaceGroupNumber().value_or(1); + const size_t matched = ApplyReferenceFreeFlags(sm.merged, static_cast(sg), + config_.reference_data); + logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)", + matched, sm.merged.size()); + } - if (result.consensus_cell && write_files && !config_.model_path.empty()) { - phase("Validating against model"); - const auto data_sg = experiment_.GetSpaceGroupNumber(); - // With a reference MTZ the merohedral indexing was already resolved against it (rotation - // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. - ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, - config_.output_prefix, logger, - data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, - /*probe_indexing_ambiguity=*/config_.reference_data.empty()); + // Radiation-damage report (rotation): the per-image scale CC-to-merge and mosaicity across the + // 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()) { + const auto &outs = indexer->GetIntegrationOutcome(); + const int nf = static_cast(outs.size()); + constexpr int nb = 10; + std::array cc_sum{}, mos_sum{}; + std::array cc_n{}, mos_n{}; + for (int f = 0; f < nf; ++f) { + const int b = nf > 1 ? std::min(nb - 1, f * nb / nf) : 0; + if (outs[f].image_scale_cc && std::isfinite(*outs[f].image_scale_cc)) { + cc_sum[b] += *outs[f].image_scale_cc; cc_n[b]++; + } + if (outs[f].mosaicity_deg && std::isfinite(*outs[f].mosaicity_deg)) { + mos_sum[b] += *outs[f].mosaicity_deg; mos_n[b]++; + } + } + int tot = 0; + for (int n : cc_n) tot += n; + if (tot >= 20) { + std::ostringstream os; + os << fmt::format("Radiation-damage report (per-image scale over {} frames):\n", nf); + os << " dose \n"; + for (int b = 0; b < nb; ++b) { + if (cc_n[b] == 0 && mos_n[b] == 0) continue; + os << fmt::format(" {:3d}-{:3d}% {:>7} {:>7}\n", b * 10, (b + 1) * 10, + cc_n[b] ? fmt::format("{:.3f}", cc_sum[b] / cc_n[b]) : std::string("-"), + mos_n[b] ? fmt::format("{:.3f}", mos_sum[b] / mos_n[b]) : std::string("-")); + } + if (cc_n[0] > 0 && cc_n[nb - 1] > 0) { + const double cc0 = cc_sum[0] / cc_n[0], cc9 = cc_sum[nb - 1] / cc_n[nb - 1]; + os << fmt::format(" => per-image CC to merge {:.3f} (first 10%) -> {:.3f} (last 10%){}", + cc0, cc9, + cc9 < cc0 - 0.05 ? " (falling: radiation damage / crystal decay)" : ""); + } + // Relative B-factor across the sweep (measured before any decay correction, against the + // low-dose start): the resolution-dependent complement to the CC/mosaicity read above. + // Radiation damage fades the high-resolution intensity, so only a POSITIVE change is dose. + // A curve no straight line describes gets no first->last number at all, and a batch whose + // data could not measure one prints "-" rather than a value. + if (!sm.statistics.radiation_damage_b_batch.empty()) { + if (std::isfinite(sm.statistics.radiation_damage_delta_b)) + os << fmt::format("\n => relative B-factor change over run = {:+.2f} A^2 (first->last){}", + sm.statistics.radiation_damage_delta_b, + sm.statistics.radiation_damage_delta_b > 5.0 + ? " (significant dose-dependent scaling: radiation damage)" : ""); + else + os << "\n => relative B-factor: the per-batch curve below is not a trend, so no " + "first->last number describes it - dose does not come back, so whatever moved " + "here was not dose; see the sweep-quality report"; + os << fmt::format("\n per-batch relative-B (A^2, {:.0f} deg/batch):", + sm.statistics.radiation_damage_batch_deg); + for (float bb : sm.statistics.radiation_damage_b_batch) + os << (std::isfinite(bb) ? fmt::format(" {:.1f}", bb) : std::string(" -")); + } + logger.Info("{}", os.str()); + result.radiation_damage_text = os.str(); + } + + // Sweep-quality report: the stretches of the sweep over which the crystal delivered much + // less than the rest of the run, and what each one looks like. Nothing is excluded because + // of it - the frames still carry signal, and this is a message for the beamline, not a + // filter. Frame numbers are processed-image ordinals, inclusive, as in _image.dat. + const auto &sq = sm.statistics.sweep_quality; + if (sq.measured) { + std::ostringstream os; + os << fmt::format("Sweep quality over {:.0f} deg (per-image scale with the incident flux, " + "which varied {:.2f}x, already divided out):\n", + sq.sweep_deg, sq.flux_peak_to_trough); + if (sq.ranges.empty()) { + os << " no stretch of the sweep is materially worse than the run"; + } else { + os << " frames rotation diagnosis severity scale CC indexed\n"; + for (const auto &r : sq.ranges) + os << fmt::format(" {:<17s} {:6.1f} deg {:<16s} {:5.2f} {:5.2f} {:5.2f} {:4.0f}%\n", + fmt::format("{}-{}", r.first_image, r.last_image), r.rotation_deg, + SweepQualityReasonText(r.reason), r.severity, r.mean_relative_scale, + r.mean_relative_cc, 100.0 * r.indexed_fraction); + os << " => severity is the fraction of the run's typical diffracting power missing " + "over the range"; + if (sq.modulation_peak_to_trough >= 1.05f) + os << fmt::format("\n => once-per-revolution modulation of the per-image scale: " + "{:.1f}x peak to trough", sq.modulation_peak_to_trough); + } + logger.Info("{}", os.str()); + } + } + + if (result.consensus_cell && write_files && config_.write_merged) { + phase("Writing reflections"); + const ErrorModelReport em_report{ + result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", + result.error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", result.error_model_isa_asymptotic) + : std::string(), + result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(), + result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()}; + WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, + em_report, result.twinning, config_.output_prefix); + // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS + // comparison. The offline self-scaling result is otherwise not exposed (process.h5's + // per-image arrays are only filled on the online per-image path). Sourced from the + // partials, which carry the first-pass per-image scale. + ScalingResult(indexer->GetIntegrationOutcome()).SaveToFile(config_.output_prefix); + } + + if (result.consensus_cell && write_files && !config_.model_path.empty()) { + phase("Validating against model"); + const auto data_sg = experiment_.GetSpaceGroupNumber(); + // With a reference MTZ the merohedral indexing was already resolved against it (rotation + // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. + ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, + config_.output_prefix, logger, + data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, + /*probe_indexing_ambiguity=*/config_.reference_data.empty()); + } } } @@ -2773,8 +2821,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.written_master_path = config_.output_prefix + "_process.h5"; } - if (observer) - observer->OnPhase(cancelled_ ? "Cancelled" : "Done"); + MarkPhase(observer, cancelled_ ? "Cancelled" : "Done"); // Total wall time for the whole run, including scaling/merging/writing above (not just the // per-image pass), so the reported rate reflects everything rugnux did. @@ -2792,5 +2839,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.used_beam_x_pxl = experiment_.GetBeamX_pxl(); result.used_beam_y_pxl = experiment_.GetBeamY_pxl(); result.used_distance_mm = experiment_.GetDetectorDistance_mm(); + + // Close the phase still running and hand over the breakdown. The timings accumulate in the object + // across passes, so the last pass's result carries the whole run rather than its own share of it. + ClosePhase(); + result.phase_timings = phase_timings_; return result; } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index bf7fd54c..dcf44b1a 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -121,6 +122,13 @@ struct ProcessConfig { std::string model_path; }; +// One phase of a run: how long it took and how much of the machine it used while it ran. +struct PhaseTiming { + std::string name; + double wall_s = 0.0; + double cpu_s = 0.0; +}; + struct ProcessResult { bool cancelled = false; uint64_t images_processed = 0; @@ -162,6 +170,11 @@ struct ProcessResult { // Twinning analysis of the final merged intensities (l_test_pairs == 0 when not computed). TwinningAnalysisResult twinning; + // Whether the space group above was reached by PROMOTING the point group during the search rather + // than being given. Carried here rather than read off `twinning`, which is only filled once the + // final merge has run - and the rotation pre-pass stops before that. + bool promoted_point_group = false; + // Space group used for the final re-scale/merge and written to the master file: determined by the // search when the user did not fix one, otherwise the fixed group. std::optional space_group_number; @@ -198,6 +211,16 @@ struct ProcessResult { // Conditions that need a person's attention, one plain sentence each (an ambiguous space group, a // symmetry axis on the spindle, an indexing ambiguity). The same messages the log warns about. std::vector warnings; + + // Wall time in each phase of the run, in the order the phases ran, summed over every pass. cpu_s is + // the CPU time the whole process burned during that phase, so cpu_s/wall_s is the mean number of + // cores the phase actually kept busy - which is what says whether a slow phase is working hard or + // sitting on one thread. Zero where the platform does not report process CPU time. + std::vector phase_timings; + + // Workers the per-image loop ran, which is NOT the thread count on a machine with GPUs: the loop is + // capped per card. The per-stage means below are worker means, so this is what they divide by. + int loop_workers = 1; }; // Callbacks for progress and live results. Methods may be called from worker threads, so an @@ -217,6 +240,16 @@ class Rugnux { ProcessConfig config_; std::atomic cancelled_{false}; + // Phase timing, accumulated across every pass of a run. MarkPhase closes the phase that was running + // and opens the named one; a phase entered twice (the two passes of a rotation run) accumulates into + // one entry rather than appearing twice, so the report reads as "where did the run go". + std::vector phase_timings_; + std::chrono::steady_clock::time_point phase_start_{}; + double phase_start_cpu_s_ = 0.0; + std::string phase_current_; + void MarkPhase(RugnuxObserver *observer, const std::string &phase); + void ClosePhase(); + // Per-frame (by image number) mosaicity fitted by RotationScaleMerge in the two-pass geometry pre-pass, // fed to the second pass's Bragg prediction. Empty outside the rotation two-pass. std::vector prepass_mosaicity_; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index efa2989a..ea464cba 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2180,7 +2180,10 @@ static int RunRugnux(int argc, char **argv) { // Dividing is a lower bound (a worker that is idle rather than blocked is not counted), so the // remainder is shown against the loop's own wall time rather than hidden. const auto &t = result.mean_processing_time; - const double per_worker = std::max(1, nthreads); + // Divide by the workers the LOOP ran, not by -N. With GPUs present the loop is capped at a few + // workers per card, so on a 2-GPU box at -N 48 those are 8 and 48 - dividing by the thread count + // understated every stage by 6x, on the one line people tune against. + const double per_worker = std::max(1, result.loop_workers); auto stage = [&](const char *name, float mean_s) { // A stage that never ran has no mean at all - the per-image indexing and scaling timers are // never fed on the two-pass rotation path, where the lattice is forced and the merge happens @@ -2188,7 +2191,7 @@ static int RunRugnux(int argc, char **argv) { return std::isfinite(mean_s) ? fmt::format(" {} {:.2f}", name, mean_s * 1e3 / per_worker) : std::string(); }; - std::cout << fmt::format("Per-image cost (ms, {} workers):", nthreads) + std::cout << fmt::format("Per-image cost (ms, {} loop workers):", result.loop_workers) << stage("decompress", t.compression) << stage("preprocess", t.preprocessing) << stage("azint", t.azint) << stage("spot-finding", t.spot_finding) << stage("indexing", t.indexing) << stage("refinement", t.refinement) @@ -2208,6 +2211,32 @@ static int RunRugnux(int argc, char **argv) { loop_ms, result.image_loop_time_s, std::max(0.0, outside_s)) << std::endl; } + // Where the whole run went, phase by phase, summed over every pass. cores = CPU time burned during + // the phase divided by its wall time, i.e. how much of the machine it actually kept busy: a phase + // with a large share and ~1 core is one thread doing all the work, which is a different problem + // from a phase that is simply large. Printed for every run so a slow dataset can be diagnosed from + // its log alone, without profiling it. + if (!result.phase_timings.empty()) { + double total = 0.0; + for (const auto &ph : result.phase_timings) total += ph.wall_s; + const bool have_cpu = std::any_of(result.phase_timings.begin(), result.phase_timings.end(), + [](const PhaseTiming &ph) { return ph.cpu_s > 0.0; }); + std::cout << "\nTime by phase (whole run, all passes):" << std::endl; + auto row = [&](const std::string &name, double wall, double cpu) { + std::string line = fmt::format(" {:<34} {:8.2f} s {:5.1f}%", name, wall, + total > 0.0 ? 100.0 * wall / total : 0.0); + if (have_cpu) + line += fmt::format(" {:6.1f} cores", wall > 0.0 ? cpu / wall : 0.0); + return line; + }; + for (const auto &ph : result.phase_timings) + if (ph.wall_s >= 0.05) + std::cout << row(ph.name, ph.wall_s, ph.cpu_s) << std::endl; + double cpu_total = 0.0; + for (const auto &ph : result.phase_timings) cpu_total += ph.cpu_s; + std::cout << row("TOTAL", total, cpu_total) << std::endl; + } + if (result.cancelled) logger.Warning("Processing was cancelled after {} images", result.images_processed);