diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index a8163887..0b979636 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -496,6 +496,12 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, }; ApplyPredictionRange(settings_prediction, experiment, latt); + // Online is bounded by what the image-buffer slot can carry, offline by the prediction limit; both + // come from BraggPrediction so the cap, the prediction and the transport headroom cannot drift. + // The predictor applies it, so what is over the cap is never integrated. + prediction.output_limit = real_time ? BraggPrediction::kOnlineMaxReflections + : BraggPrediction::kPredictionOutput; + // Predict, then integrate with the selected integrator (box-sum or profile-fit). auto pred_start_time = std::chrono::steady_clock::now(); auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction); @@ -511,25 +517,6 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, auto integration_end_time = std::chrono::steady_clock::now(); msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); - // Online is bounded by what the image-buffer slot can carry; offline by the prediction limit. Both - // come from BraggPrediction so the cap, the prediction and the transport headroom cannot drift. - const size_t kMaxReflections = real_time ? BraggPrediction::kOnlineMaxReflections - : BraggPrediction::kPredictionOutput; - if (i_outcome.reflections.size() > kMaxReflections) { - // Keep only smallest d (highest resolution) - std::nth_element(i_outcome.reflections.begin(), - i_outcome.reflections.begin() + static_cast(kMaxReflections), - i_outcome.reflections.end(), - [](const Reflection& a, const Reflection& b) { - return a.d < b.d; - }); - i_outcome.reflections.resize(kMaxReflections); - - // Optional: make output ordered by d (nice for downstream / debugging) - std::sort(i_outcome.reflections.begin(), i_outcome.reflections.end(), - [](const Reflection& a, const Reflection& b) { return a.d < b.d; }); - } - CalcISigma(msg, i_outcome.reflections); CalcWilsonBFactor(msg, i_outcome.reflections); diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index f18a95be..87e042bc 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -14,18 +14,24 @@ void BraggPrediction::GrowCapacity(int count) { } int BraggPrediction::TruncateToOutput(int count) { - if (count <= kPredictionOutput) + if (count <= output_limit) return count; - std::partial_sort(reflections.begin(), reflections.begin() + kPredictionOutput, + // Partiality first, then excitation error: exactly one of the two says anything on each path. On the + // rotation path dist_ewald is identically zero - BraggPredictionRot picks the rocking coordinate so + // that 2*S0.p + p.p = 0, i.e. |S| = 1/lambda exactly - so ranking by it there left an hkl-lexicographic + // prefix rather than the best-recorded reflections. On the still path every partiality is 1, so the + // ranking falls through to dist_ewald, which is what it has always been. + std::partial_sort(reflections.begin(), reflections.begin() + output_limit, reflections.begin() + count, [](const Reflection &a, const Reflection &b) { + if (a.partiality != b.partiality) return a.partiality > b.partiality; if (a.dist_ewald != b.dist_ewald) return a.dist_ewald < b.dist_ewald; if (a.h != b.h) return a.h < b.h; if (a.k != b.k) return a.k < b.k; return a.l < b.l; }); - return kPredictionOutput; + return output_limit; } namespace { diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index e9b9dfb6..1c634ce3 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -47,13 +47,14 @@ protected: // different merged reflections on a CPU-only build than on a GPU one. virtual void GrowCapacity(int count); - // Deterministically cap Calc's output at kPredictionOutput: if more were predicted, keep the ones - // closest to the Ewald sphere (smallest excitation error), ties broken by hkl. Returns the kept - // count. Below the cap it is a no-op. Call at the end of every Calc override. + // Deterministically cap Calc's output at output_limit: if more were predicted, keep the best-recorded + // ones - largest partiality on the rotation path, and, since partiality is 1 for every still, smallest + // excitation error on the still path - with hkl breaking what is left. Returns the kept count. Below + // the cap it is a no-op. Call at the end of every Calc override. int TruncateToOutput(int count); public: // The prediction buffer holds up to kPredictionCapacity reflections so a strong lattice does not - // overflow it. Calc returns at most kPredictionOutput, the number that flows downstream and is serialized - kept low so the + // overflow it. Calc returns at most output_limit, the number that flows downstream and is serialized - kept low so the // per-image reflection list stays within the frame transport headroom. // Starting size. On the GPU path the buffer grows to whatever a frame actually predicts // (GrowCapacity), so a large cell is not truncated there; the CPU path still caps at this value, @@ -61,12 +62,12 @@ public: // and NON-DETERMINISTIC - the GPU kernels claim slots with an atomicAdd, so which reflections // survived depended on block scheduling and changed between runs of the same command. static constexpr int kPredictionCapacity = 20000; - // How many reflections may flow downstream per image. Sized for a large unit cell: a ~2.8e6 A^3 - // cell predicts up to ~44000 per frame at 2.4 A. Truncating below what the frame really has costs - // more than it saves - the selection keeps the smallest excitation errors, i.e. the nearly - // fully-recorded reflections, and the rotation combine rebuilds a full FROM the partials it drops - // (measured on such a crystal: CC1/2 98 -> 60, ISa 8.6 -> 1.8). DiffractionExperiment's image-buffer - // headroom is derived from this, so the transport can carry what the analysis produces. + // How many reflections may flow downstream per image, offline. Sized for a large unit cell: a + // ~2.8e6 A^3 cell predicts up to ~44000 per frame at 2.4 A. Truncating below what the frame really + // has costs more than it saves - the selection keeps the best-recorded reflections, and the rotation + // combine rebuilds a full FROM the partials it drops (measured on such a crystal: CC1/2 98 -> 60, + // ISa 8.6 -> 1.8). DiffractionExperiment's image-buffer headroom is derived from this, so the + // transport can carry what the analysis produces. static constexpr int kPredictionOutput = 65536; // What the ONLINE path may carry per image. The acquisition transports every reflection list through // a fixed-size image-buffer slot, and the slot size divides a fixed total - so sizing the slot for @@ -75,6 +76,11 @@ public: // has no such budget, keeps the full one. DiffractionExperiment's buffer headroom derives from THIS. static constexpr int kOnlineMaxReflections = 10000; + // How many reflections Calc may return. A caller with a tighter cap than the offline one - online, + // whose transport slot holds kOnlineMaxReflections - sets its own here, so the surplus is dropped + // before it is integrated rather than integrated and then thrown away. + int output_limit = kPredictionOutput; + explicit BraggPrediction(int max_reflections = kPredictionCapacity); virtual ~BraggPrediction() = default; diff --git a/image_analysis/indexing/IndexerThreadPool.cpp b/image_analysis/indexing/IndexerThreadPool.cpp index ef59725e..ee5d11e3 100644 --- a/image_analysis/indexing/IndexerThreadPool.cpp +++ b/image_analysis/indexing/IndexerThreadPool.cpp @@ -14,10 +14,28 @@ #include "FFTIndexerCPU.h" #endif -IndexerThread::IndexerThread(const IndexingSettings &settings, int threadid) { +// The indexer for one RESOLVED algorithm, or nullptr if this build/host cannot serve it. +static std::unique_ptr MakeIndexer(IndexingAlgorithmEnum algorithm, const IndexingSettings &settings) { +#ifdef JFJOCH_USE_CUDA + if (get_gpu_count() > 0) { + if (algorithm == IndexingAlgorithmEnum::FFT) + return std::make_unique(settings); + if (algorithm == IndexingAlgorithmEnum::FFBIDX) + return std::make_unique(); + } +#endif +#ifdef JFJOCH_USE_FFTW + if (algorithm == IndexingAlgorithmEnum::FFTW) + return std::make_unique(settings); +#endif + return nullptr; +} + +IndexerThread::IndexerThread(const IndexingSettings &settings, int threadid, IndexerConstruction construction) + : settings_(settings), construction_(construction) { std::unique_lock lock(m); state = TaskState::STARTING; - worker_thread = std::thread(&IndexerThread::Worker, this, std::cref(settings), threadid); + worker_thread = std::thread(&IndexerThread::Worker, this, threadid); c_running.wait(lock, [this] { return state != TaskState::STARTING; }); if (state == TaskState::ERROR) { worker_thread.join(); @@ -26,7 +44,7 @@ IndexerThread::IndexerThread(const IndexingSettings &settings, int threadid) { } } -void IndexerThread::Worker(const IndexingSettings &settings, int threadid) { +void IndexerThread::Worker(int threadid) { try { pin_gpu(); } catch (const std::exception &e) { @@ -37,58 +55,42 @@ void IndexerThread::Worker(const IndexingSettings &settings, int threadid) { std::unique_ptr fft_indexer, ffbidx_indexer, fftw_indexer; -#ifdef JFJOCH_USE_CUDA - try { - if (get_gpu_count() > 0) { - if (settings.GetAlgorithm() == IndexingAlgorithmEnum::Auto - || settings.GetAlgorithm() == IndexingAlgorithmEnum::FFT) - fft_indexer = std::make_unique(settings); + // Preconstruct: build every indexer the requested algorithm could resolve to before the pool + // reports ready, so no cuFFT planning happens once frames are flowing, and a failure is fatal + // for the pool instead of being met frame by frame. OnFirstUse skips this and builds in the + // dispatch below. + if (construction_ == IndexerConstruction::Preconstruct) { + try { + const auto requested = settings_.GetAlgorithm(); - if (settings.GetAlgorithm() == IndexingAlgorithmEnum::Auto - || settings.GetAlgorithm() == IndexingAlgorithmEnum::FFBIDX) - ffbidx_indexer = std::make_unique(); + if (requested == IndexingAlgorithmEnum::Auto || requested == IndexingAlgorithmEnum::FFT) + fft_indexer = MakeIndexer(IndexingAlgorithmEnum::FFT, settings_); + + if (requested == IndexingAlgorithmEnum::Auto || requested == IndexingAlgorithmEnum::FFBIDX) + ffbidx_indexer = MakeIndexer(IndexingAlgorithmEnum::FFBIDX, settings_); + + if ((requested == IndexingAlgorithmEnum::Auto && get_gpu_count() == 0) + || requested == IndexingAlgorithmEnum::FFTW) + fftw_indexer = MakeIndexer(IndexingAlgorithmEnum::FFTW, settings_); + } catch (const std::exception &e) { + spdlog::error("Failed to initialize indexer: {}", e.what()); + { + std::unique_lock lock(m); + state = TaskState::ERROR; + } + c_running.notify_all(); + return; + } catch (...) { + spdlog::error("Failed to initialize indexer"); + { + std::unique_lock lock(m); + state = TaskState::ERROR; + } + c_running.notify_all(); + return; } - } catch (const std::exception &e) { - spdlog::error("Failed to initialize GPU indexer: {}", e.what()); - { - std::unique_lock lock(m); - state = TaskState::ERROR; - } - c_running.notify_all(); - return; - } catch (...) { - spdlog::error("Failed to initialize GPU indexer"); - { - std::unique_lock lock(m); - state = TaskState::ERROR; - } - c_running.notify_all(); - return; } -#endif -#ifdef JFJOCH_USE_FFTW - try { - if ((settings.GetAlgorithm() == IndexingAlgorithmEnum::Auto && (get_gpu_count() == 0)) - || settings.GetAlgorithm() == IndexingAlgorithmEnum::FFTW) - fftw_indexer = std::make_unique(settings); - } catch (const std::exception &e) { - spdlog::error("Failed to initialize FFTW indexer: {}", e.what()); - { - std::unique_lock lock(m); - state = TaskState::ERROR; - } - c_running.notify_all(); - return; - } catch (...) { - spdlog::error("Failed to initialize FFTW indexer"); - { - std::unique_lock lock(m); - state = TaskState::ERROR; - } - c_running.notify_all(); - return; - } -#endif + { std::unique_lock lock(m); state = TaskState::IDLE; @@ -110,28 +112,33 @@ void IndexerThread::Worker(const IndexingSettings &settings, int threadid) { std::unique_ptr tmp_result; try { auto algorithm = input->experiment.GetIndexingAlgorithm(); - Indexer *indexer = nullptr; + std::unique_ptr *slot = nullptr; - if (algorithm == IndexingAlgorithmEnum::FFT && fft_indexer) { - indexer = fft_indexer.get(); - } else if (algorithm == IndexingAlgorithmEnum::FFBIDX && ffbidx_indexer) { - indexer = ffbidx_indexer.get(); - } else if (algorithm == IndexingAlgorithmEnum::FFTW && fftw_indexer) { - indexer = fftw_indexer.get(); - } else { + switch (algorithm) { + case IndexingAlgorithmEnum::FFT: slot = &fft_indexer; break; + case IndexingAlgorithmEnum::FFBIDX: slot = &ffbidx_indexer; break; + case IndexingAlgorithmEnum::FFTW: slot = &fftw_indexer; break; + default: break; + } + + // A preconstructing worker already holds it; an OnFirstUse worker builds it here, + // on the first frame that resolves to this algorithm. + if (slot && !*slot) + *slot = MakeIndexer(algorithm, settings_); + + if (!slot || !*slot) { // Algorithm is already resolved here (never Auto/None - see - // IndexerThreadPool::Run). Reaching this means the resolved algorithm - // has no matching indexer in this worker (e.g. a GPU algorithm on a - // host without a GPU) - fail loudly instead of silently not indexing. + // IndexerThreadPool::Run, which also checked this host can serve it). Reaching + // this means the resolved algorithm has no matching indexer in this build - + // fail loudly instead of silently not indexing. throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Internal error: no indexer available for the resolved " "indexing algorithm"); } - if (indexer) { - indexer->Setup(input->experiment); - tmp_result = std::make_unique(indexer->Run(input->recip)); - } + Indexer &indexer = **slot; + indexer.Setup(input->experiment); + tmp_result = std::make_unique(indexer.Run(input->recip)); } catch (std::exception &e) { tmp_result = nullptr; spdlog::error("Indexer thread {} failed: {}", threadid, e.what()); @@ -182,13 +189,13 @@ IndexerThread::~IndexerThread() { Finalize(); } -IndexerThreadPool::IndexerThreadPool(const IndexingSettings &settings) +IndexerThreadPool::IndexerThreadPool(const IndexingSettings &settings, IndexerConstruction construction) : worker_busy(settings.GetIndexingThreads(), 0), worker_free_count(settings.GetIndexingThreads()), viable_cell_min_spots(settings.GetViableCellMinSpots()), blocking(settings.GetBlockingBehavior()) { for (size_t i = 0; i < settings.GetIndexingThreads(); ++i) - tasks.emplace_back(std::make_unique(std::cref(settings), i)); + tasks.emplace_back(std::make_unique(std::cref(settings), i, construction)); } int IndexerThreadPool::GetFreeWorker() { diff --git a/image_analysis/indexing/IndexerThreadPool.h b/image_analysis/indexing/IndexerThreadPool.h index fd7e0a04..8d980683 100644 --- a/image_analysis/indexing/IndexerThreadPool.h +++ b/image_analysis/indexing/IndexerThreadPool.h @@ -19,12 +19,34 @@ #include "../common/DiffractionExperiment.h" #include "Indexer.h" +// When a worker builds its indexers. +// +// Preconstruct (the default, and what the online service needs): every indexer the requested +// algorithm could resolve to is built while the pool is being constructed, so it is resident and +// planned before the first frame arrives. Building one means a cuFFT plan plus ~0.6 GB of device +// allocation, and the first cuFFT call in a process also pays the library's one-time init (0.3 s +// measured); jfjoch_broker cannot have any of that land on the live data path, so holding an indexer +// the resolved algorithm may never dispatch to is an accepted cost there. +// +// OnFirstUse (offline batch - rugnux, jfjoch_viewer): nothing is latency-critical, so the indexer is +// built on the first frame that needs it. The algorithm is only resolved from the frame's +// DiffractionExperiment (rotation -> FFT, stills with a known cell -> FFBIDX), so preconstructing +// leaves each worker holding a fully allocated FFTIndexerGPU - 0.6 GB of cuFFT plan and histograms - +// that can never be dispatched to. +enum class IndexerConstruction { Preconstruct, OnFirstUse }; + class IndexerThread { struct TaskInput { const DiffractionExperiment &experiment; const std::vector &recip; }; + // Held by value: with IndexerConstruction::OnFirstUse the worker builds its indexer long after + // the constructor returned, and pools are routinely built from a temporary - for instance + // IndexerThreadPool(experiment.GetIndexingSettings()), which returns by value. + const IndexingSettings settings_; + const IndexerConstruction construction_; + bool stop = false; enum class TaskState {STARTING, IDLE, READY, RUNNING, COMPLETED, ERROR} state = TaskState::STARTING; std::mutex m; @@ -35,9 +57,9 @@ class IndexerThread { std::unique_ptr task_input = nullptr; std::thread worker_thread; - void Worker(const IndexingSettings& settings, int threadid); + void Worker(int threadid); public: - IndexerThread(const IndexingSettings& settings, int threadid); + IndexerThread(const IndexingSettings& settings, int threadid, IndexerConstruction construction); ~IndexerThread(); std::unique_ptr Run(const DiffractionExperiment &experiment, const std::vector &recip); void Finalize(); @@ -53,7 +75,8 @@ class IndexerThreadPool { const bool blocking; int GetFreeWorker(); public: - IndexerThreadPool(const IndexingSettings& settings); + IndexerThreadPool(const IndexingSettings& settings, + IndexerConstruction construction = IndexerConstruction::Preconstruct); IndexerResult Run(const DiffractionExperiment& experiment, const std::vector& recip); }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index e509dcda..afae34ac 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -109,7 +109,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ // First-pass engines, built from the current (nominal) geometry: an index-only pass (no // integration) to obtain each frame's spots + assigned HKL + orientation. AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); - IndexerThreadPool pool(experiment_.GetIndexingSettings()); + IndexerThreadPool pool(experiment_.GetIndexingSettings(), IndexerConstruction::OnFirstUse); IndexAndRefine indexer(experiment_, &pool, /*retain_outcomes=*/false); auto pass_settings = config_.spot_finding; @@ -523,7 +523,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // never uses. Stills fire the FFT once per image across all workers, so keep the default there. if (config_.rotation_indexing && config_.two_pass_rotation && experiment_.IsRotationIndexing()) indexing_settings.IndexingThreads(2); - indexer_pool = std::make_unique(indexing_settings); + indexer_pool = std::make_unique(indexing_settings, IndexerConstruction::OnFirstUse); indexer = std::make_unique(experiment_, indexer_pool.get()); if (!config_.reference_data.empty()) indexer->ReferenceIntensities(config_.reference_data); @@ -762,7 +762,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b auto coarse_settings = experiment_.GetIndexingSettings(); coarse_settings.FFT_HighResolution_A(3.5f); // low-order reflections only -> robust long axis coarse_settings.IndexingThreads(2); - IndexerThreadPool coarse_pool(coarse_settings); + IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse); // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full- // resolution per-frame validation would be low (the coarse cell is metrically right but // imprecise), so do NOT score it here - it is only the reference for the constrained re-index. @@ -799,9 +799,29 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } if (!cancelled_) { + // The two ways the first pass can end without a usable lattice. They are one situation from + // the user's side - no crystal was found in these frames - so both name the cause and then + // give the same remedies. + constexpr const char *next_steps = + " Nothing was integrated. Check the beam centre (--beam-x / --beam-y), raise " + "--max-spots if the frames are well populated with spots, and use --azint-only to " + "write the per-image analysis for inspection."; if (!best.result.has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Two-pass rotation indexing failed"); + std::string("Two-pass rotation indexing found no lattice.") + next_steps); + // A lattice that indexes almost none of the validation frames is not this crystal's lattice. + // Forcing it anyway integrates every image against it; the few partials that survive yield no + // fulls in the 3D combine, and the run then dies inside the merge complaining about resolution - + // far from the cause. Stop here instead, while the reason is still visible. The floor is one + // frame in six, several times below the weakest crystal that still merges. It is checked AFTER + // the long-axis rescue, so a metric that rescue recovers is never rejected on its pre-rescue + // score. + if (best.score < static_cast(validation.size()) / 6) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + fmt::format("Two-pass rotation indexing found only a lattice that " + "indexes {}/{} validation frames - too few for it to be " + "this crystal's lattice.{}", + best.score, static_cast(validation.size()), next_steps)); indexer->ForceRotationIndexerResult(*best.result); logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames)", best.name, best.score, static_cast(validation.size())); diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 51907f29..94a5a5da 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -106,7 +106,7 @@ JFJochImageReadingWorker::JFJochImageReadingWorker(const SpotFindingSettings &se qRegisterMetaType("ReferenceMtzInfo"); spot_finding_settings = settings; - indexing = std::make_unique(indexing_settings); + indexing = std::make_unique(indexing_settings, IndexerConstruction::OnFirstUse); http_reader.Experiment(experiment); file_reader.Experiment(experiment); thumb_color_scale_.Select(ColorScaleEnum::Indigo);