diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d557b1749..a1401382a 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -871,6 +871,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b saturation = std::async(std::launch::async, [&] { return SampleSaturation(start_image, images_to_process); }); if (full && config_.refine_geometry.has_value() && !config_.detector_distance_unknown) RefineStillsGeometry(start_image, end_image, images_to_process, observer); + // The detector tilt's fits read only the frames just fitted and the geometry they left: started now, they + // run beside the pre-scan and the geometry trial, and are used where neither moved the geometry. + if (full) + StartStillsTiltAhead(); // The geometry pre-pass integrates only the first part of the sweep (pass 2 reads it all), and on a // rotation sweep the decisions before pass 2 read only that part too (PrepassEnd). @@ -890,13 +894,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // The stills geometry refinement above ran from the file's centre and distance, within a few pixels and // millimetres of them; where the frames have just moved either, it is run again from where they are. + bool refitted = false; if (full && !experiment_.IsRotationIndexing() && config_.beam_center_check && !config_.live - && DecideStillsGeometry(start_image, images_to_process) && config_.refine_geometry.has_value()) + && DecideStillsGeometry(start_image, images_to_process) && config_.refine_geometry.has_value()) { RefineStillsGeometry(start_image, end_image, images_to_process, observer); + refitted = true; + } // The detector tilt, now that the centre is settled: a centre error on frames indexed at it is exactly a // tilt's pattern shift, so a tilt fitted before the centre is decided takes up the centre's error. if (full) - DecideStillsTilt(); + DecideStillsTilt(refitted); // The geometry pre-pass integrates at the radius the run started with; the width the pre-scan just // measured is applied to the canonical pass instead (RunAllPasses, between the two). See diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index d99a804c3..20d350226 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -1418,8 +1418,13 @@ class Rugnux { int64_t on_lattice = 0; // ... of which on the lattice }; std::vector FindStillsSampleSpots(const std::vector &ordinals, int start_image); - StillsTrial IndexStillsSample(const DiffractionExperiment &x, IndexerThreadPool &pool, - const std::vector &frames); + // A geometry hypothesis and the frames it is indexed on. + struct StillsHypothesis { + DiffractionExperiment x; + const std::vector *frames; + }; + // Every hypothesis indexed on its frames, all frames of all hypotheses side by side. + std::vector IndexStillsSample(const std::vector &hypotheses, IndexerThreadPool &pool); // Returns whether the run's geometry was moved off the file's. bool DecideStillsGeometry(int start_image, int images_to_process); // Stills: the detector tilt as a hypothesis, refined or held, kept where the frames it was not fitted to @@ -1431,7 +1436,24 @@ class Rugnux { std::vector modules; }; std::optional stills_tilt_bundle_; - void DecideStillsTilt(); + // The tilt hypothesis tried on a bundle: the held-out residuals with the tilt held and refined, the tilted fit + // on every frame, and the geometry and cell the fits started from. + struct StillsTiltTrial { + DiffractionGeometry start; + UnitCell cell; + gemmi::CrystalSystem crystal_system; + double held = 0, refined = 0; + GeometryRefinerResult tilted; + }; + // The same fits started as soon as the first RefineStillsGeometry made the bundle, beside the pre-scan and the + // geometry trial; DecideStillsTilt takes them where nothing in between moved the geometry or the bundle. + static StillsTiltTrial TryStillsTilt(const StillsTiltBundle &bundle, const DiffractionGeometry &start, + const UnitCell &cell, GeometryRefinerSettings st); + GeometryRefinerSettings StillsTiltSettings() const; + std::shared_future stills_tilt_ahead_; + void StartStillsTiltAhead(); + // bundle_refitted: RefineStillsGeometry ran again after StartStillsTiltAhead. + void DecideStillsTilt(bool bundle_refitted); // Stills: the cell the user gave (-C / reference), kept apart from the one the geometry refinement // writes into experiment_. See RefineStillsGeometry. std::optional stills_user_cell_; diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index eb5962d08..79394851f 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -115,6 +115,13 @@ GeomRefineFrame ToRefineFrame(const DataMessage &msg, int ordinal) { return f; } +bool SameGeometry(const DiffractionGeometry &a, const DiffractionGeometry &b) { + return a.GetBeamX_pxl() == b.GetBeamX_pxl() && a.GetBeamY_pxl() == b.GetBeamY_pxl() + && a.GetDetectorDistance_mm() == b.GetDetectorDistance_mm() && a.GetPixelSize_mm() == b.GetPixelSize_mm() + && a.GetWavelength_A() == b.GetWavelength_A() && a.GetPoniRot1_rad() == b.GetPoniRot1_rad() + && a.GetPoniRot2_rad() == b.GetPoniRot2_rad() && a.GetPoniRot3_rad() == b.GetPoniRot3_rad(); +} + // Spots piled up at one radius - an ice or other powder ring - are not this crystal's lattice. They are // found by their radius in millimetres about (beam_x, beam_y), not by a d-spacing: a d computed at a // distance that may be wrong misses the rings, and then they are most of the spots a frame has. @@ -554,23 +561,10 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ // back to where it was given, and refitted with the tilt on every frame they gave a worse merge than held. Where the tilt is real it bends the pattern in a way nothing // else in the fit can (the keystone); where it is not, the gain on frames the fit has not seen is nothing and the // tilt found is too small to matter either way, so no threshold has to say how large a tilt must be. -void Rugnux::DecideStillsTilt() { - if (!stills_tilt_bundle_) - return; - Logger logger("Rugnux"); - TimingMark mark(logger, "stills tilt decision"); - const StillsTiltBundle bundle = std::move(*stills_tilt_bundle_); - stills_tilt_bundle_.reset(); - - const DiffractionGeometry start = experiment_.GetDiffractionGeometry(); - const UnitCell cell = experiment_.GetUnitCell().value(); - GeometryRefinerSettings st; - st.crystal_system = experiment_.GetCrystalSystem(); - st.num_threads = config_.nthreads; - st.refine_distance = !config_.detector_distance_given; - st.cell_reg_coeff = 0.0; +// The fits read only what they are given, so they can run beside the steps before the decision (StartStillsTiltAhead). +Rugnux::StillsTiltTrial Rugnux::TryStillsTilt(const StillsTiltBundle &bundle, const DiffractionGeometry &start, + const UnitCell &cell, GeometryRefinerSettings st) { st.modules = bundle.modules; - std::vector half[2]; for (size_t i = 0; i < bundle.frames.size(); i++) half[i % 2].push_back(bundle.frames[i]); @@ -604,14 +598,58 @@ void Rugnux::DecideStillsTilt() { sf.refine_tilt = true; return RefineGlobalGeometry(start, cell, bundle.frames, sf); }); - const double held = (held_fits[0].get() + held_fits[1].get()) / 2.0; - const double refined = (tilt_fits[0].get() + tilt_fits[1].get()) / 2.0; - const GeometryRefinerResult r = tilted.get(); - const bool keep = r.ok && refined < held; + StillsTiltTrial t{start, cell, st.crystal_system}; + t.held = (held_fits[0].get() + held_fits[1].get()) / 2.0; + t.refined = (tilt_fits[0].get() + tilt_fits[1].get()) / 2.0; + t.tilted = tilted.get(); + return t; +} + +GeometryRefinerSettings Rugnux::StillsTiltSettings() const { + GeometryRefinerSettings st; + st.crystal_system = experiment_.GetCrystalSystem(); + st.num_threads = config_.nthreads; + st.refine_distance = !config_.detector_distance_given; + st.cell_reg_coeff = 0.0; + return st; +} + +void Rugnux::StartStillsTiltAhead() { + if (!stills_tilt_bundle_) + return; + stills_tilt_ahead_ = std::async(std::launch::async, TryStillsTilt, *stills_tilt_bundle_, + experiment_.GetDiffractionGeometry(), experiment_.GetUnitCell().value(), + StillsTiltSettings()).share(); +} + +void Rugnux::DecideStillsTilt(bool bundle_refitted) { + if (!stills_tilt_bundle_) + return; + Logger logger("Rugnux"); + TimingMark mark(logger, "stills tilt decision"); + std::optional ahead; + if (stills_tilt_ahead_.valid()) { + ahead = stills_tilt_ahead_.get(); + stills_tilt_ahead_ = {}; + } + const StillsTiltBundle bundle = std::move(*stills_tilt_bundle_); + stills_tilt_bundle_.reset(); + + const DiffractionGeometry start = experiment_.GetDiffractionGeometry(); + const UnitCell cell = experiment_.GetUnitCell().value(); + const GeometryRefinerSettings st = StillsTiltSettings(); + + // The fits started ahead are these fits where they started from the same bundle, geometry and cell. + const bool use_ahead = ahead && !bundle_refitted && SameGeometry(ahead->start, start) && ahead->cell == cell + && ahead->crystal_system == st.crystal_system; + logger.Debug("Detector tilt: the fits started beside the pre-scan are {}", use_ahead ? "used" : "not used"); + const StillsTiltTrial t = use_ahead ? std::move(*ahead) : TryStillsTilt(bundle, start, cell, st); + const GeometryRefinerResult &r = t.tilted; + const bool keep = r.ok && t.refined < t.held; logger.Info("Detector tilt: rot1 {:+.3f} deg, rot2 {:+.3f} deg refined ({:+.3f}, {:+.3f} held); held-out spot " "residual (capped rms) {:.3f} px held, {:.3f} px refined, {} frames - the tilt is {}", r.rot1_rad * 180.0 / PI, r.rot2_rad * 180.0 / PI, start.GetPoniRot1_rad() * 180.0 / PI, - start.GetPoniRot2_rad() * 180.0 / PI, held, refined, bundle.frames.size(), + start.GetPoniRot2_rad() * 180.0 / PI, t.held, t.refined, bundle.frames.size(), keep ? "refined" : "held"); if (!keep) return; @@ -670,46 +708,57 @@ std::vector Rugnux::FindStillsSampleSpots(const std::vector &o return frames; } -Rugnux::StillsTrial Rugnux::IndexStillsSample(const DiffractionExperiment &x, IndexerThreadPool &pool, - const std::vector &frames) { - IndexAndRefine indexer(x, &pool, /*retain_outcomes=*/false); - const DiffractionGeometry geom = x.GetDiffractionGeometry(); +std::vector Rugnux::IndexStillsSample(const std::vector &hypotheses, + IndexerThreadPool &pool) { auto settings = config_.spot_finding; settings.enable = true; settings.indexing = true; + std::vector> indexers; + std::vector> items; // (hypothesis, frame) + for (size_t h = 0; h < hypotheses.size(); h++) { + indexers.push_back(std::make_unique(hypotheses[h].x, &pool, /*retain_outcomes=*/false)); + for (size_t i = 0; i < hypotheses[h].frames->size(); i++) + items.emplace_back(h, i); + } + std::mutex m; - StillsTrial total; + std::vector total(hypotheses.size()); std::atomic next = 0; auto worker = [&]() { - StillsTrial t; - for (size_t i = next.fetch_add(1); i < frames.size() && !cancelled_; i = next.fetch_add(1)) { + std::vector t(hypotheses.size()); + for (size_t n = next.fetch_add(1); n < items.size() && !cancelled_; n = next.fetch_add(1)) { + const auto [h, i] = items[n]; + const DiffractionGeometry &geom = hypotheses[h].x.GetDiffractionGeometry(); // The spots where the frame put them, at this trial's geometry. - DataMessage msg = frames[i]; + DataMessage msg = (*hypotheses[h].frames)[i]; for (auto &s : msg.spots) { s.d_A = geom.PxlToRes(s.x, s.y); s.indexed = false; s.lattice = -1; } - t.frames++; - if (!indexer.IndexFrameOnly(msg, settings)) + t[h].frames++; + if (!indexers[h]->IndexFrameOnly(msg, settings)) continue; - t.indexed++; - t.spots += static_cast(msg.spots.size()); - t.on_lattice += msg.spot_count_indexed.value_or(0); + t[h].indexed++; + t[h].spots += static_cast(msg.spots.size()); + t[h].on_lattice += msg.spot_count_indexed.value_or(0); } - indexer.RethrowIndexerError(); std::lock_guard lock(m); - total.frames += t.frames; - total.indexed += t.indexed; - total.spots += t.spots; - total.on_lattice += t.on_lattice; + for (size_t h = 0; h < t.size(); h++) { + total[h].frames += t[h].frames; + total[h].indexed += t[h].indexed; + total[h].spots += t[h].spots; + total[h].on_lattice += t[h].on_lattice; + } }; std::vector> futures; for (int i = 0; i < config_.nthreads; ++i) futures.push_back(std::async(std::launch::async, worker)); for (auto &f : futures) f.get(); + for (const auto &indexer : indexers) + indexer->RethrowIndexerError(); return total; } @@ -785,34 +834,61 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { scan_frames.assign(frames.begin(), frames.begin() + std::min(frames.size(), scan_count)); }; find_sample(); - IndexerThreadPool pool(experiment_.GetIndexingSettings(), IndexerConstruction::OnFirstUse, gpu_work_beside_); + // As many indexer threads as the image pass has (RunPipeline): the hypotheses are indexed side by side, and + // with the default four the workers waited on the pool. + IndexingSettings pool_settings = experiment_.GetIndexingSettings(); + pool_settings.IndexingThreads(std::clamp(config_.nthreads, 1, 8)); + IndexerThreadPool pool(pool_settings, IndexerConstruction::OnFirstUse, gpu_work_beside_); // A hypothesis is indexed once: the searches below meet the same one more than once - the decided centre - // is the file's whenever the file's is kept - and it gives the same count every time. - std::map *>, int64_t> tried; - const auto trial = [&](float x, float y, float distance, const std::vector &on) { - const auto key = std::make_tuple(x, y, distance, &on); - if (const auto it = tried.find(key); it != tried.end()) - return it->second; - DiffractionExperiment e = experiment_; - place(e, x, y, distance); - if (stills_user_cell_) - e.SetUnitCell(*stills_user_cell_); // the grid asks the user's cell, not the one refined at the file's distance - // The distance scan's frames are indexed with a wider acceptance: the scan only has to find the peak, and - // it looks at centres a few pixels off and at distances up to half a step off, either of which leaves the - // strict tolerance with nothing on a large cell at a long distance (a 92 A cell at 199 mm indexed on no - // centre 3 px off). Everything after the scan is counted at the run's own tolerance. - if (&on == &scan_frames) { - IndexingSettings wide = e.GetIndexingSettings(); - wide.Tolerance(SCAN_TOLERANCE); - e.ImportIndexingSettings(wide); + // is the file's whenever the file's is kept - and it gives the same count every time. The hypotheses a + // search asks are indexed together (try_all) - a hypothesis' count does not depend on what is indexed + // beside it - and then read one by one, in the order the search compares them. + using TrialKey = std::tuple *>; + std::map tried; + // Hypotheses indexed in one go: enough to keep every worker and the indexer pool busy on 8 frames each, + // few enough that their experiments and indexers stay small. + constexpr size_t HYPOTHESES_AT_ONCE = 64; + const auto try_all = [&](const std::vector &keys) { + std::vector todo; + for (const auto &key : keys) + if (!tried.contains(key) && std::find(todo.begin(), todo.end(), key) == todo.end()) + todo.push_back(key); + for (size_t c0 = 0; c0 < todo.size(); c0 += HYPOTHESES_AT_ONCE) { + const size_t c1 = std::min(todo.size(), c0 + HYPOTHESES_AT_ONCE); + std::vector hypotheses; + for (size_t k = c0; k < c1; k++) { + const auto &[x, y, distance, on] = todo[k]; + DiffractionExperiment e = experiment_; + place(e, x, y, distance); + if (stills_user_cell_) + e.SetUnitCell(*stills_user_cell_); // the grid asks the user's cell, not the one refined at the file's distance + // The distance scan's frames are indexed with a wider acceptance: the scan only has to find the peak, and + // it looks at centres a few pixels off and at distances up to half a step off, either of which leaves the + // strict tolerance with nothing on a large cell at a long distance (a 92 A cell at 199 mm indexed on no + // centre 3 px off). Everything after the scan is counted at the run's own tolerance. + if (on == &scan_frames) { + IndexingSettings wide = e.GetIndexingSettings(); + wide.Tolerance(SCAN_TOLERANCE); + e.ImportIndexingSettings(wide); + } + hypotheses.push_back(StillsHypothesis{std::move(e), on}); + } + const std::vector results = IndexStillsSample(hypotheses, pool); + for (size_t k = c0; k < c1; k++) { + const auto &[x, y, distance, on] = todo[k]; + const StillsTrial &t = results[k - c0]; + logger.Debug("Stills geometry trial: centre ({:.2f},{:.2f}) distance {:.2f} mm - {}/{} frames indexed, " + "{} of their {} spots on the lattice", x, y, distance, t.indexed, t.frames, + t.on_lattice, t.spots); + tried[todo[k]] = t.on_lattice; + } } - const StillsTrial t = IndexStillsSample(e, pool, on); - logger.Debug("Stills geometry trial: centre ({:.2f},{:.2f}) distance {:.2f} mm - {}/{} frames indexed, " - "{} of their {} spots on the lattice", x, y, distance, t.indexed, t.frames, - t.on_lattice, t.spots); - tried[key] = t.on_lattice; - return t.on_lattice; + }; + const auto trial = [&](float x, float y, float distance, const std::vector &on) { + const TrialKey key{x, y, distance, &on}; + try_all({key}); + return tried.at(key); }; // The centres: the file's, the measured one and the point-symmetry capture's - two estimators of the @@ -826,6 +902,10 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { float x = file_x, y = file_y; int64_t score_here = 0; const auto decide_centre = [&](float distance) { + std::vector keys{{file_x, file_y, distance, &frames}}; + for (const auto &[cx, cy] : candidates) + keys.emplace_back(cx, cy, distance, &frames); + try_all(keys); const int64_t file_score = trial(file_x, file_y, distance, frames); int64_t best = file_score; float best_x = file_x, best_y = file_y; @@ -867,6 +947,10 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { float best_distance = file_distance, best_x = x, best_y = y; int64_t best = score_here; const auto fine = [&](float about, float cx, float cy) { + std::vector keys; + for (const float f : {0.98f, 0.99f, 1.0f, 1.01f, 1.02f}) + keys.emplace_back(cx, cy, about * f, &frames); + try_all(keys); for (const float f : {0.98f, 0.99f, 1.0f, 1.01f, 1.02f}) { const int64_t score = trial(cx, cy, about * f, frames); if (score > best) { @@ -893,6 +977,11 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { centres.emplace_back(file_x, file_y); // The local maxima of every centre's scan, strongest first; a distance is looked at once. std::vector> maxima; + std::vector keys; + for (const auto &[cx, cy] : centres) + for (const float d : grid) + keys.emplace_back(cx, cy, d, &scan_frames); + try_all(keys); for (const auto &[cx, cy] : centres) { std::vector scan(grid.size()); for (size_t i = 0; i < grid.size(); i++) @@ -909,6 +998,12 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { float grid_best = file_distance, grid_x = x, grid_y = y; int64_t grid_score = -1; + keys.clear(); + for (const size_t p : peaks) + for (size_t i = (p > 0 ? p - 1 : 0); i <= std::min(grid.size() - 1, p + 1); i++) + for (const auto &[cx, cy] : centres) + keys.emplace_back(cx, cy, grid[i], &strongest); + try_all(keys); for (const size_t p : peaks) for (size_t i = (p > 0 ? p - 1 : 0); i <= std::min(grid.size() - 1, p + 1); i++) for (const auto &[cx, cy] : centres) { @@ -977,6 +1072,12 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { int64_t g_best = trial(x, y, distance, strongest); for (const auto &[step, reach] : {std::pair{3.0f, 4}, std::pair{1.0f, 2}}) { const float cx0 = gx, cy0 = gy; + std::vector keys; + for (int i = -reach; i <= reach; i++) + for (int j = -reach; j <= reach; j++) + if (i != 0 || j != 0) + keys.emplace_back(cx0 + step * i, cy0 + step * j, distance, &strongest); + try_all(keys); for (int i = -reach; i <= reach; i++) for (int j = -reach; j <= reach; j++) { if (i == 0 && j == 0) @@ -989,6 +1090,7 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { } } } + try_all({{x, y, distance, &frames}, {gx, gy, distance, &frames}}); const int64_t here = trial(x, y, distance, frames), there = trial(gx, gy, distance, frames); logger.Info("Stills beam centre grid at {:.2f} mm: ({:.2f},{:.2f}) puts {} sampled spots on the lattice, " "({:.2f},{:.2f}) {}", distance, gx, gy, there, x, y, here);