Merge fc-trial: the stills geometry trial indexes its hypotheses side by side; the tilt fits run beside the pre-scan (decision-identical)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu

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

rugnux/Rugnux.h

@@ -1418,8 +1418,13 @@ class Rugnux {
         int64_t on_lattice = 0;    // ... of which on the lattice
     };
     std::vector<DataMessage> FindStillsSampleSpots(const std::vector<int> &ordinals, int start_image);
-    StillsTrial IndexStillsSample(const DiffractionExperiment &x, IndexerThreadPool &pool,
-                                  const std::vector<DataMessage> &frames);
+    // A geometry hypothesis and the frames it is indexed on.
+    struct StillsHypothesis {
+        DiffractionExperiment x;
+        const std::vector<DataMessage> *frames;
+    };
+    // Every hypothesis indexed on its frames, all frames of all hypotheses side by side.
+    std::vector<StillsTrial> IndexStillsSample(const std::vector<StillsHypothesis> &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<ModuleShift> modules;
     };
     std::optional<StillsTiltBundle> 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<StillsTiltTrial> 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<UnitCell> stills_user_cell_;

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<GeomRefineFrame> 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<StillsTiltTrial> 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<DataMessage> Rugnux::FindStillsSampleSpots(const std::vector<int> &o
     return frames;
 }
 
-Rugnux::StillsTrial Rugnux::IndexStillsSample(const DiffractionExperiment &x, IndexerThreadPool &pool,
-                                              const std::vector<DataMessage> &frames) {
-    IndexAndRefine indexer(x, &pool, /*retain_outcomes=*/false);
-    const DiffractionGeometry geom = x.GetDiffractionGeometry();
+std::vector<Rugnux::StillsTrial> Rugnux::IndexStillsSample(const std::vector<StillsHypothesis> &hypotheses,
+                                                           IndexerThreadPool &pool) {
     auto settings = config_.spot_finding;
     settings.enable = true;
     settings.indexing = true;
 
+    std::vector<std::unique_ptr<IndexAndRefine>> indexers;
+    std::vector<std::pair<size_t, size_t>> items;   // (hypothesis, frame)
+    for (size_t h = 0; h < hypotheses.size(); h++) {
+        indexers.push_back(std::make_unique<IndexAndRefine>(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<StillsTrial> total(hypotheses.size());
     std::atomic<size_t> 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<StillsTrial> 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<int64_t>(msg.spots.size());
-            t.on_lattice += msg.spot_count_indexed.value_or(0);
+            t[h].indexed++;
+            t[h].spots += static_cast<int64_t>(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<std::future<void>> 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<size_t>(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<std::tuple<float, float, float, const std::vector<DataMessage> *>, int64_t> tried;
-    const auto trial = [&](float x, float y, float distance, const std::vector<DataMessage> &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<float, float, float, const std::vector<DataMessage> *>;
+    std::map<TrialKey, int64_t> 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<TrialKey> &keys) {
+        std::vector<TrialKey> 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<StillsHypothesis> 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<StillsTrial> 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<DataMessage> &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<TrialKey> 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<TrialKey> 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<std::pair<int64_t, size_t>> maxima;
+            std::vector<TrialKey> 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<int64_t> 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<TrialKey> 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);