diff --git a/docs/JFJOCH_BROKER.md b/docs/JFJOCH_BROKER.md index 7ea6040e2..7fb6bde35 100644 --- a/docs/JFJOCH_BROKER.md +++ b/docs/JFJOCH_BROKER.md @@ -266,15 +266,11 @@ all: on fewer, too few index for the geometry refinement to run, and the stated buffer has to hold them, and the frames that arrive while the pre-pass runs, in half its slots (for a JUNGFRAU 8M, a slot is the size of an uncompressed frame); a smaller buffer starts the pre-pass on as many as it holds. -Even 500 frames give the pre-pass only a few dozen indexed ones at a serial hit rate, and none where the -first frames are empty (a fixed-target chip whose crystals sit further in). So, where a cell is given and the -pre-pass refined the geometry on fewer than 200 frames, the centre and distance are refined again on the frames -as they are analysed, from their spots and lattices, without holding their images: once 200 strong frames (20 -indexed spots or more) are in, or 5000 indexed, the geometry is fitted, as rugnux fits it on its sample, and the -analysis switches to it; a second round does the same on frames analysed at that geometry. Two fits in all, a -fixed cost. The frames analysed before the last switch are written as analysed but left out of the merge (their -reflections were predicted at the earlier geometry, and the images are no longer there to integrate again); -the switches are notes of the live status, and the end message carries the geometry the run ended at. +The geometry is fixed when the pre-pass ends and does not change during the collection, so every frame is analysed +and merged at it and the radial profile has one binning throughout. Where the pre-pass could not refine it (too few +of its frames indexed, as with a small image buffer, a low hit rate or a chip whose first frames are empty, or no unit +cell given), the distance in the header is used for the whole run: this is a warning in the log and a note of the +live status, and a wrong header distance then costs resolution and indexed frames. Frames arriving during the pre-pass, or faster than the analysis, wait in the image buffer, in at most half of its slots; a frame that finds no free slot is dropped and counted as for a full image buffer, so analysis diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 2bbfcb9ad..7dc673a16 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -1128,12 +1128,6 @@ IntegrationOutcome IndexAndRefine::TakeIntegrationOutcome(int64_t image) { return std::move(integration_outcome.at(image)); } -void IndexAndRefine::DropIntegrationOutcome(int64_t image) { - const std::unique_lock ul(reflections_mutex); - integration_outcome.at(image) = {}; - unit_cells.at(image).reset(); -} - std::vector &IndexAndRefine::GetIntegrationOutcome() { return integration_outcome; } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 0467e31ad..f90c3267b 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -309,7 +309,5 @@ public: // The retained outcome of one image, moved out: what is left in its place is empty. Safe beside the // analysis of other images. IntegrationOutcome TakeIntegrationOutcome(int64_t image); - // The retained outcome of one image dropped, and its cell with it: it is not in the consensus. - void DropIntegrationOutcome(int64_t image); const std::vector &GetIntegrationOutcome() const; }; diff --git a/receiver/StillsGPUAnalysis.cpp b/receiver/StillsGPUAnalysis.cpp index 9eac35441..5ebeda6f3 100644 --- a/receiver/StillsGPUAnalysis.cpp +++ b/receiver/StillsGPUAnalysis.cpp @@ -207,15 +207,15 @@ void StillsGPUAnalysis::PrePassThread() { setup.mapping = std::make_unique(setup.experiment, setup.mask); // Retained for the merge, which takes each image's reflections as soon as it is analysed. indexer = std::make_unique(setup.experiment, indexing_thread_pool, merge_frames, true); - // The pre-pass's fit rests on the frames of its sample that indexed; where those are fewer than the - // running refinement bundles, the geometry is refined again, from where the pre-pass left it, on the - // frames as they are analysed. - if (geometry_frames < StillsRunningGeometry::STRONG_FRAMES && setup.experiment.GetUnitCell()) { - running_geometry = std::make_unique(setup.experiment, MERGE_THREADS); - frame_epoch.assign(experiment.GetImageNum(), 0); - geometry_state = GEOMETRY_PENDING; - logger.Info("The pre-pass refined the geometry on {} frames: it is refined again on {} strong frames as " - "they are analysed", geometry_frames, StillsRunningGeometry::STRONG_FRAMES); + if (geometry_frames == 0) { + const std::string text = fmt::format( + "The pre-pass could not refine the geometry (too few frames indexed, or no unit cell given): the " + "whole run is analysed at the distance in the header, {:.2f} mm, beam centre ({:.2f}, {:.2f}) px", + setup.experiment.GetDetectorDistance_mm(), setup.experiment.GetBeamX_pxl(), + setup.experiment.GetBeamY_pxl()); + logger.Warning("{}", text); + if (status) + status->Note(text); } if (merge_frames) { merge_experiment = setup.experiment; @@ -267,54 +267,32 @@ void StillsGPUAnalysis::AnalysisThread() { return; } - // A thread whose analysis cannot be built still passes its frames on, unanalysed. It is built again at - // the geometry the running refinement switches to. + // A thread whose analysis cannot be built still passes its frames on, unanalysed. std::unique_ptr analysis; - int epoch = -1; - const auto build = [&] { - epoch = geometry_epoch; - analysis.reset(); - try { - analysis = std::make_unique(setup.experiment, *setup.mapping, setup.mask, - *indexer, /*enable_fused_adaptive_gpu=*/true); - } catch (const std::exception &e) { - AnalysisError(std::string("analysis thread could not start, its frames are written unanalysed - ") - + e.what()); - } - }; + try { + analysis = std::make_unique(setup.experiment, *setup.mapping, setup.mask, + *indexer, /*enable_fused_adaptive_gpu=*/true); + } catch (const std::exception &e) { + AnalysisError(std::string("analysis thread could not start, its frames are written unanalysed - ") + + e.what()); + } while (auto held = queue.GetBlocking()) { DataMessage &message = held->message; std::optional profile; - bool switch_geometry = false; const auto start_time = std::chrono::steady_clock::now(); - { - std::unique_lock gate(gate_mutex); - gate_cv.wait(gate, [this] { return !switch_waiting; }); - } - { - std::shared_lock sl(geometry_mutex); - if (epoch != geometry_epoch) - build(); - profile.emplace(*setup.mapping); - try { - if (analysis) - analysis->Analyze(message, *profile, settings()); - } catch (const JFJochException &e) { - // The frame is still written, without its analysis. - logger.ErrorException(e); - } catch (const std::exception &e) { - AnalysisError(fmt::format("image {} written unanalysed - {}", message.number, e.what())); - } - if (geometry_state == GEOMETRY_PENDING) { - frame_epoch.at(message.number) = static_cast(epoch + 1); - switch_geometry = running_geometry->Add(message); - } + profile.emplace(*setup.mapping); + try { + if (analysis) + analysis->Analyze(message, *profile, settings()); + } catch (const JFJochException &e) { + // The frame is still written, without its analysis. + logger.ErrorException(e); + } catch (const std::exception &e) { + AnalysisError(fmt::format("image {} written unanalysed - {}", message.number, e.what())); } message.processing_time_s = std::chrono::duration(std::chrono::steady_clock::now() - start_time).count(); - if (switch_geometry) - SwitchGeometry(); try { sink(message, *profile); if (merge) @@ -337,63 +315,6 @@ void StillsGPUAnalysis::AnalysisThread() { } } -void MoveStillsGeometry(DiffractionExperiment &experiment, float beam_x_pxl, float beam_y_pxl, float distance_mm) { - AzimuthalIntegrationSettings azint = experiment.GetAzimuthalIntegrationSettings(); - azint.QRange_recipA(azint.GetLowQ_recipA(), azint.GetHighQ_recipA()); - experiment.ImportAzimuthalIntegrationSettings(azint); - experiment.BeamX_pxl(beam_x_pxl).BeamY_pxl(beam_y_pxl).DetectorDistance_mm(distance_mm); -} - -void StillsGPUAnalysis::SwitchGeometry() { - // One fit, in the analysis thread that offered the last frame it needed; the others go on meanwhile. - const auto t0 = std::chrono::steady_clock::now(); - const auto r = running_geometry->Refine(); - const bool last_round = ++geometry_round == GEOMETRY_ROUNDS; - if (!r) { - logger.Warning("Running geometry refinement, round {}: no fit on {} frames ({} strong); the frames are " - "analysed at the geometry they were", geometry_round, running_geometry->Offered(), - running_geometry->Strong()); - geometry_state = geometry_epoch > 0 ? GEOMETRY_SWITCHED : GEOMETRY_KEPT; - return; - } - const std::string text = fmt::format( - "Geometry refined on {} analysed frames ({} strong), round {}: beam centre ({:.2f}, {:.2f}) -> " - "({:.2f}, {:.2f}) px, distance {:.2f} -> {:.2f} mm ({:.1f} s)", running_geometry->Offered(), - running_geometry->Strong(), geometry_round, setup.experiment.GetBeamX_pxl(), - setup.experiment.GetBeamY_pxl(), r->beam_x_px, r->beam_y_px, setup.experiment.GetDetectorDistance_mm(), - r->distance_mm, std::chrono::duration(std::chrono::steady_clock::now() - t0).count()); - { - std::unique_lock gate(gate_mutex); - switch_waiting = true; - } - { - std::unique_lock ul(geometry_mutex); - MoveStillsGeometry(setup.experiment, static_cast(r->beam_x_px), static_cast(r->beam_y_px), - static_cast(r->distance_mm)); - retired_mappings.push_back(std::move(setup.mapping)); - setup.mapping = std::make_unique(setup.experiment, setup.mask); - ++geometry_epoch; - if (last_round) - geometry_state = GEOMETRY_SWITCHED; - else - running_geometry = std::make_unique(setup.experiment, MERGE_THREADS); - } - { - std::unique_lock gate(gate_mutex); - switch_waiting = false; - } - gate_cv.notify_all(); - { - std::unique_lock ul(live_mutex); - live.beam_x_pxl = r->beam_x_px; - live.beam_y_pxl = r->beam_y_px; - live.distance_mm = r->distance_mm; - } - logger.Info("{}", text); - if (status) - status->Note(text); -} - std::optional StillsGPUAnalysis::GetConsensusUnitCell() const { return indexer->GetConsensusUnitCell(); } @@ -446,27 +367,8 @@ void StillsGPUAnalysis::MergeThread() { } else if (next > 0 && merge->Added() % PREPARE_EVERY == 0) merge->Prepare(MERGE_THREADS); }; - // The frames analysed at a geometry the running refinement may still replace wait until it is decided: - // left out where analysed before its last switch, merged as analysed otherwise. - std::vector held_back; - const auto release = [&] { - for (const int64_t h: held_back) - if (frame_epoch[h] - 1 < geometry_epoch) { - indexer->DropIntegrationOutcome(h); - add(h, IntegrationOutcome{}); - } else - add(h, indexer->TakeIntegrationOutcome(h)); - held_back.clear(); - }; - for (int64_t n = merge_queue.GetBlocking(); n >= 0; n = merge_queue.GetBlocking()) { - if (!frame_epoch.empty() && frame_epoch[n] > 0) - held_back.push_back(n); - else - add(n, indexer->TakeIntegrationOutcome(n)); - if (geometry_state != GEOMETRY_PENDING && !held_back.empty()) - release(); - } - release(); + for (int64_t n = merge_queue.GetBlocking(); n >= 0; n = merge_queue.GetBlocking()) + add(n, indexer->TakeIntegrationOutcome(n)); if (have_cell()) { merged = merge->Final(static_cast(AvailableCpus())); logger.Info("Stills merge after the last image: {:.2f} s", merged->seconds); diff --git a/receiver/StillsGPUAnalysis.h b/receiver/StillsGPUAnalysis.h index d37b8148c..071662810 100644 --- a/receiver/StillsGPUAnalysis.h +++ b/receiver/StillsGPUAnalysis.h @@ -4,7 +4,6 @@ #pragma once #include -#include #include #include #include @@ -13,7 +12,6 @@ #include #include #include -#include #include #include "StillsPrePass.h" @@ -26,7 +24,6 @@ #include "../common/ThreadSafeFIFO.h" #include "../image_analysis/IndexAndRefine.h" #include "../rugnux/StillsLiveMerge.h" -#include "../rugnux/StillsRunningGeometry.h" // Serial frames analysed on the GPU as the DECTRIS ones are (MXAnalysisWithoutFPGA), after a pre-pass on // the first PREPASS_IMAGES (RunStillsPrePass) has decided the beam centre and the mask they are analysed @@ -44,17 +41,10 @@ // before Finish returns (GetMerged). The cell is the run's consensus at the first merge, as rugnux merges in // it (the collection's cell picks the frames merged); the space group the collection's, else P1. // merge = false leaves it out. -// Where the pre-pass refined the centre and distance on fewer frames than StillsRunningGeometry bundles, or -// could not refine them (too few of the frames it had indexed: a small image buffer, a serial hit rate, a chip -// whose first frames are empty), and a cell is given, they are refined again while the frames are analysed -// (StillsRunningGeometry, from the spots and lattices of the frames analysed, no image held), and the analysis -// switches to that geometry as each of GEOMETRY_ROUNDS fits is made. The frames analysed before the last -// switch are not merged (their reflections were predicted at an earlier geometry; the images are gone); where -// no fit can be made they are merged as analysed. -// Moves the analysis to a refined geometry during the run. The q range stays where it was, also one left to -// reach the detector's edge: the start message and the run's profile are binned on it. -void MoveStillsGeometry(DiffractionExperiment &experiment, float beam_x_pxl, float beam_y_pxl, float distance_mm); - +// The geometry is fixed when the pre-pass ends: every frame is analysed and merged at it. Where the pre-pass could +// not refine it (too few of the frames it had indexed: a small image buffer, a serial hit rate, a chip whose first +// frames are empty; or no cell given), the header's distance stays for the whole run, with a warning in the log +// and the status. class StillsGPUAnalysis { public: // Enough frames for the distance to be refined at a serial hit rate: of 100 frames of a collection indexed @@ -147,29 +137,6 @@ private: std::optional merged; std::future merge_future; - // The geometry refined while the frames are analysed (see above), in GEOMETRY_ROUNDS rounds, each on the - // frames analysed at the last one's geometry (the spots of a frame are indexed at the geometry it was - // analysed at, and a start 2.6% off left the first fit 0.6% off). The analysis takes geometry_mutex shared - // for each frame; a switch takes it exclusive and raises geometry_epoch, which the analysis threads see to - // rebuild their analysis. geometry_state is PENDING until the last round, then SWITCHED or KEPT. - // frame_epoch: per frame number, 1 + the epoch it was analysed at while PENDING (0: after); the merge - // leaves out those analysed before the last switch. - constexpr static int GEOMETRY_ROUNDS = 2; - enum GeometryState : int { GEOMETRY_PENDING = 0, GEOMETRY_SWITCHED = 1, GEOMETRY_KEPT = 2 }; - std::unique_ptr running_geometry; - int geometry_round = 0; - std::shared_mutex geometry_mutex; - // A switch waiting for geometry_mutex closes this gate to new frames, or the frames that keep taking the - // mutex shared would hold it off to the end of the collection. - std::mutex gate_mutex; - std::condition_variable gate_cv; - bool switch_waiting = false; - std::atomic geometry_epoch{0}; - std::atomic geometry_state{GEOMETRY_KEPT}; - std::vector frame_epoch; - std::vector> retired_mappings; // the analyses before the switch read them - void SwitchGeometry(); - void StartPrePass(); void PrePassThread(); void AnalysisThread(); diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 72ab81c3d..bd1462509 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -25,8 +25,6 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxLive.h StillsLiveMerge.cpp StillsLiveMerge.h - StillsRunningGeometry.cpp - StillsRunningGeometry.h RugnuxRefine.cpp RugnuxRefine.h Rugnux.h diff --git a/rugnux/StillsRunningGeometry.cpp b/rugnux/StillsRunningGeometry.cpp deleted file mode 100644 index e5be725dc..000000000 --- a/rugnux/StillsRunningGeometry.cpp +++ /dev/null @@ -1,95 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "StillsRunningGeometry.h" - -#include - -#include "../image_analysis/scale_merge/CellPopulation.h" - -namespace { - // A frame the geometry refinement can use, as rugnux takes it from its sample: the spots on lattice 0 with - // their index, the others without. - GeomRefineFrame RefineFrame(const DataMessage &msg) { - GeomRefineFrame f; - f.lattice = *msg.indexing_lattice; - f.ordinal = static_cast(msg.number); - for (const auto &s: msg.spots) { - if (s.indexed && s.lattice == 0) - f.spots.push_back(GeomRefineSpot{s.x, s.y, static_cast(s.h), static_cast(s.k), - static_cast(s.l), s.panel}); - else - f.unindexed.push_back(GeomRefineSpot{s.x, s.y, 0, 0, 0, s.panel}); - } - return f; - } -} - -StillsRunningGeometry::StillsRunningGeometry(const DiffractionExperiment &in_experiment, int in_nthreads) - : experiment(in_experiment), nthreads(in_nthreads) {} - -bool StillsRunningGeometry::Add(const DataMessage &msg) { - std::optional f; - if (msg.indexing_result.value_or(false) && msg.indexing_lattice.has_value()) - f = RefineFrame(msg); - std::unique_lock ul(m); - if (done) - return false; - ++offered; - if (f && f->spots.size() >= 6) { - ++indexed; - strong += static_cast(f->spots.size()) >= MIN_STRONG_SPOTS; - frames.push_back(std::move(*f)); - } - done = static_cast(strong) >= STRONG_FRAMES || indexed >= INDEXED_BUDGET; - return done; -} - -int64_t StillsRunningGeometry::Offered() const { - std::unique_lock ul(m); - return offered; -} - -size_t StillsRunningGeometry::Strong() const { - std::unique_lock ul(m); - return strong; -} - -std::optional StillsRunningGeometry::Refine() { - std::unique_lock ul(m); - const auto cell = experiment.GetUnitCell(); - if (!cell) - return {}; - // As RugnuxStills.cpp's sample: only the frames of the given cell's population are fitted, the population - // told from the strong frames. - const double given = CrystalLattice(*cell).CalcVolume(); - std::vector strong_cells; - for (const auto &f: frames) - if (static_cast(f.spots.size()) >= MIN_STRONG_SPOTS) - strong_cells.push_back(f.lattice.GetUnitCell()); - if (strong_cells.empty()) - return {}; - const auto [low, high] = GivenCellPopulation(strong_cells, *cell); - std::vector bundle; - for (auto &f: frames) { - const double r = f.lattice.CalcVolume() / given - 1.0; - if (static_cast(f.spots.size()) >= MIN_STRONG_SPOTS && r >= low && r <= high) - bundle.push_back(std::move(f)); - } - frames.clear(); - std::sort(bundle.begin(), bundle.end(), [](const GeomRefineFrame &a, const GeomRefineFrame &b) { - if (a.spots.size() != b.spots.size()) return a.spots.size() > b.spots.size(); - return a.ordinal < b.ordinal; - }); - if (static_cast(bundle.size()) > STRONG_FRAMES) - bundle.resize(STRONG_FRAMES); - - GeometryRefinerSettings settings; - settings.crystal_system = experiment.GetCrystalSystem(); - settings.num_threads = nthreads; - settings.cell_reg_coeff = 0.0; // cell and distance both free, as rugnux fits them - const GeometryRefinerResult r = RefineGlobalGeometry(experiment.GetDiffractionGeometry(), *cell, bundle, settings); - if (!r.ok) - return {}; - return r; -} diff --git a/rugnux/StillsRunningGeometry.h b/rugnux/StillsRunningGeometry.h deleted file mode 100644 index ad3260673..000000000 --- a/rugnux/StillsRunningGeometry.h +++ /dev/null @@ -1,47 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include - -#include "../common/DiffractionExperiment.h" -#include "../common/JFJochMessages.h" -#include "../image_analysis/geom_refinement/GeometryRefiner.h" - -// The detector centre and distance of a serial collection refined from frames as they are analysed, holding -// no image: the live path's geometry where its pre-pass could not refine it (too few of the frames the image -// buffer holds indexed). The frames' indexed spots and lattices are kept until STRONG_FRAMES strong frames -// (at least MIN_STRONG_SPOTS indexed spots) are in, or INDEXED_BUDGET frames indexed (empty frames - a chip's -// wells without a crystal - do not count); then the strongest -// are bundled and the geometry refined once, as rugnux refines it on its sample (cell and distance free, the -// detector untilted), on the frames of the given cell's population. A fixed cost: one bundle, one fit. -// Thread-safe. -class StillsRunningGeometry { -public: - constexpr static int STRONG_FRAMES = 200; - constexpr static int MIN_STRONG_SPOTS = 20; - constexpr static int64_t INDEXED_BUDGET = 5000; - - // `experiment`: the geometry the frames are analysed at, the cell given. - explicit StillsRunningGeometry(const DiffractionExperiment &experiment, int nthreads); - - // A frame as analysed. True once: for the frame after which Refine is to be called. - bool Add(const DataMessage &msg); - // The refined geometry, where the fit ran. - std::optional Refine(); - // Frames offered and frames kept, for the log - [[nodiscard]] int64_t Offered() const; - [[nodiscard]] size_t Strong() const; - -private: - const DiffractionExperiment experiment; - const int nthreads; - mutable std::mutex m; - std::vector frames; - int64_t offered = 0, indexed = 0; - size_t strong = 0; - bool done = false; -}; diff --git a/tests/StillsGPUAnalysisTest.cpp b/tests/StillsGPUAnalysisTest.cpp index b890ac222..ef803a2e7 100644 --- a/tests/StillsGPUAnalysisTest.cpp +++ b/tests/StillsGPUAnalysisTest.cpp @@ -125,30 +125,3 @@ TEST_CASE("StillsGPUAnalysis_FailingFrameDoesNotStopAnalysis", "[StillsGPUAnalys CHECK(s.steps.back().duration_s); CHECK(buffer.GetStatus().preparation_slots == 0); } - -TEST_CASE("StillsGPUAnalysis_MovedGeometryKeepsTheProfileBins", "[StillsGPUAnalysis]") { - DiffractionExperiment x = StillsExperiment(); - const PixelMask mask(x); - const AzimuthalIntegrationMapping before(x, mask); - AzimuthalIntegrationProfile run_profile(before); - - // Left to itself the high q follows the detector's edge, and with it the number of bins. - DiffractionExperiment unfrozen = x; - unfrozen.DetectorDistance_mm(120); - REQUIRE(AzimuthalIntegrationMapping(unfrozen, mask).GetBinNumber() != before.GetBinNumber()); - - MoveStillsGeometry(x, 510, 260, 120); - CHECK(x.GetDetectorDistance_mm() == Catch::Approx(120)); - CHECK(x.GetBeamX_pxl() == Catch::Approx(510)); - CHECK(x.GetBeamY_pxl() == Catch::Approx(260)); - - const AzimuthalIntegrationMapping after(x, mask); - CHECK(after.GetBinNumber() == before.GetBinNumber()); - CHECK(after.GetBinToQ() == before.GetBinToQ()); - const AzimuthalIntegrationProfile frame_profile(after); - CHECK_NOTHROW(run_profile += frame_profile); - - // A second move stays on the same bins. - MoveStillsGeometry(x, 505, 255, 90); - CHECK(AzimuthalIntegrationMapping(x, mask).GetBinToQ() == before.GetBinToQ()); -}