indexing: select predicted reflections by partiality, build indexers where it pays
Build Packages / build:viewer-tgz:cpu (push) Successful in 7m34s
Build Packages / build:viewer-tgz:cuda (push) Successful in 8m42s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 13m24s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 13m31s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 13m44s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 14m5s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 14m16s
Build Packages / build:rpm (rocky8) (push) Successful in 11m28s
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 12m45s
Build Packages / XDS test (durin plugin) (push) Successful in 7m39s
Build Packages / Generate python client (push) Successful in 36s
Build Packages / Build documentation (push) Successful in 1m4s
Build Packages / Create release (push) Skipped
Build Packages / build:rpm (rocky9) (push) Successful in 12m20s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 12m35s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 13m9s
Build Packages / DIALS test (push) Successful in 13m57s
Build Packages / XDS test (neggia plugin) (push) Successful in 7m57s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 8m39s
Build Packages / Unit tests (push) Successful in 1h1m5s
Build Packages / build:windows:nocuda (push) Failing after 2s
Build Packages / build:windows:cuda (push) Failing after 3s
Build Packages / build:viewer-tgz:cpu (push) Successful in 7m34s
Build Packages / build:viewer-tgz:cuda (push) Successful in 8m42s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 13m24s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 13m31s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 13m44s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 14m5s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 14m16s
Build Packages / build:rpm (rocky8) (push) Successful in 11m28s
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 12m45s
Build Packages / XDS test (durin plugin) (push) Successful in 7m39s
Build Packages / Generate python client (push) Successful in 36s
Build Packages / Build documentation (push) Successful in 1m4s
Build Packages / Create release (push) Skipped
Build Packages / build:rpm (rocky9) (push) Successful in 12m20s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 12m35s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 13m9s
Build Packages / DIALS test (push) Successful in 13m57s
Build Packages / XDS test (neggia plugin) (push) Successful in 7m57s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 8m39s
Build Packages / Unit tests (push) Successful in 1h1m5s
Build Packages / build:windows:nocuda (push) Failing after 2s
Build Packages / build:windows:cuda (push) Failing after 3s
When more reflections are predicted for a frame than the output can hold, the surplus was dropped by keeping those closest to the Ewald sphere. On the rotation path that quantity is identically zero by construction - the rocking coordinate is chosen so the scattering vector lands exactly on the sphere - so the comparison fell through to h, k and l and the survivors were whichever came first in lexicographic order. Measured on a large cell: every value within one float ulp of zero, and the kept set had a MEAN PARTIALITY BELOW that of the full set, i.e. worse than choosing at random. Rank by partiality instead, which the predictor already computes and which is what the header always claimed was being kept. On the one regression crystal large enough to cross the cap this lifts completeness from 84.8% to 90.2% on the same observations; multiplicity and R_meas move the way they must when the same measurements cover more of reciprocal space. The online path asked for a cap of ten thousand but the truncation was hardcoded to the offline limit, so the broker predicted and integrated up to six times what it could transport and discarded the rest after paying for it. Honour the caller's limit, which also makes the post-integration re-truncation dead code. Indexer pool construction becomes a policy. The online service needs every indexer resident before data arrives, because a cuFFT plan built on the first frame is planning time inside the measurement; spending memory to be ready is the intended trade there and stays the default. Offline there is no such deadline, and a stills run with a known cell was holding a fully allocated FFT indexer per worker that the algorithm resolution can never dispatch - 2.8 GB where 0.4 GB is needed. rugnux and the viewer opt into building on first use; the broker, the receiver and the tests are untouched. This also removes a dangling reference that was latent: the worker held the settings by reference although the pool is routinely constructed from a temporary, which only survived because eager construction finished inside the constructor call. Finally, refuse a first-pass lattice that indexes fewer than a sixth of the validation frames. It fires on nothing in the regression set - the weakest real crystal sits at 22 of 60, more than twice the floor - so it is a backstop, but the failure it prevents is one the set does contain: a dataset with no crystal at all adopts a lattice from its powder rings, integrates every image against it, and dies much later inside the merge complaining about resolution. It now stops in the first pass and says what to try. Regression set: 36 of 37 crystals byte-identical, the exception being the completeness gain above; 34 of 37 space groups, no failures. Full unit suite passes. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -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<float>(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<long>(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);
|
||||
|
||||
|
||||
@@ -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 {
|
||||
|
||||
@@ -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;
|
||||
|
||||
@@ -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<Indexer> MakeIndexer(IndexingAlgorithmEnum algorithm, const IndexingSettings &settings) {
|
||||
#ifdef JFJOCH_USE_CUDA
|
||||
if (get_gpu_count() > 0) {
|
||||
if (algorithm == IndexingAlgorithmEnum::FFT)
|
||||
return std::make_unique<FFTIndexerGPU>(settings);
|
||||
if (algorithm == IndexingAlgorithmEnum::FFBIDX)
|
||||
return std::make_unique<FFBIDXIndexer>();
|
||||
}
|
||||
#endif
|
||||
#ifdef JFJOCH_USE_FFTW
|
||||
if (algorithm == IndexingAlgorithmEnum::FFTW)
|
||||
return std::make_unique<FFTIndexerCPU>(settings);
|
||||
#endif
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
IndexerThread::IndexerThread(const IndexingSettings &settings, int threadid, IndexerConstruction construction)
|
||||
: settings_(settings), construction_(construction) {
|
||||
std::unique_lock<std::mutex> 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<Indexer> 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<FFTIndexerGPU>(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<FFBIDXIndexer>();
|
||||
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<std::mutex> lock(m);
|
||||
state = TaskState::ERROR;
|
||||
}
|
||||
c_running.notify_all();
|
||||
return;
|
||||
} catch (...) {
|
||||
spdlog::error("Failed to initialize indexer");
|
||||
{
|
||||
std::unique_lock<std::mutex> 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<std::mutex> lock(m);
|
||||
state = TaskState::ERROR;
|
||||
}
|
||||
c_running.notify_all();
|
||||
return;
|
||||
} catch (...) {
|
||||
spdlog::error("Failed to initialize GPU indexer");
|
||||
{
|
||||
std::unique_lock<std::mutex> 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<FFTIndexerCPU>(settings);
|
||||
} catch (const std::exception &e) {
|
||||
spdlog::error("Failed to initialize FFTW indexer: {}", e.what());
|
||||
{
|
||||
std::unique_lock<std::mutex> lock(m);
|
||||
state = TaskState::ERROR;
|
||||
}
|
||||
c_running.notify_all();
|
||||
return;
|
||||
} catch (...) {
|
||||
spdlog::error("Failed to initialize FFTW indexer");
|
||||
{
|
||||
std::unique_lock<std::mutex> lock(m);
|
||||
state = TaskState::ERROR;
|
||||
}
|
||||
c_running.notify_all();
|
||||
return;
|
||||
}
|
||||
#endif
|
||||
|
||||
{
|
||||
std::unique_lock<std::mutex> lock(m);
|
||||
state = TaskState::IDLE;
|
||||
@@ -110,28 +112,33 @@ void IndexerThread::Worker(const IndexingSettings &settings, int threadid) {
|
||||
std::unique_ptr<IndexerResult> tmp_result;
|
||||
try {
|
||||
auto algorithm = input->experiment.GetIndexingAlgorithm();
|
||||
Indexer *indexer = nullptr;
|
||||
std::unique_ptr<Indexer> *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<IndexerResult>(indexer->Run(input->recip));
|
||||
}
|
||||
Indexer &indexer = **slot;
|
||||
indexer.Setup(input->experiment);
|
||||
tmp_result = std::make_unique<IndexerResult>(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<IndexerThread>(std::cref(settings), i));
|
||||
tasks.emplace_back(std::make_unique<IndexerThread>(std::cref(settings), i, construction));
|
||||
}
|
||||
|
||||
int IndexerThreadPool::GetFreeWorker() {
|
||||
|
||||
@@ -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<Coord> &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<TaskInput> 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<IndexerResult> Run(const DiffractionExperiment &experiment, const std::vector<Coord> &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<Coord>& recip);
|
||||
};
|
||||
|
||||
|
||||
+24
-4
@@ -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<IndexerThreadPool>(indexing_settings);
|
||||
indexer_pool = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse);
|
||||
indexer = std::make_unique<IndexAndRefine>(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<int>(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<int>(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<int>(validation.size()));
|
||||
|
||||
@@ -106,7 +106,7 @@ JFJochImageReadingWorker::JFJochImageReadingWorker(const SpotFindingSettings &se
|
||||
qRegisterMetaType<ReferenceMtzInfo>("ReferenceMtzInfo");
|
||||
spot_finding_settings = settings;
|
||||
|
||||
indexing = std::make_unique<IndexerThreadPool>(indexing_settings);
|
||||
indexing = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse);
|
||||
http_reader.Experiment(experiment);
|
||||
file_reader.Experiment(experiment);
|
||||
thumb_color_scale_.Select(ColorScaleEnum::Indigo);
|
||||
|
||||
Reference in New Issue
Block a user