Bragg prediction: one limit per index, not one cube

Each Miller index is bounded by its OWN axis - |h| <= a/d_min, |k| <= b/d_min,
|l| <= c/d_min - so a single half-width has to be sized for the longest axis and
then walks the short ones far past anything the resolution cut can keep. Give the
predictor max_h, max_k and max_l instead, in all four implementations (CPU and GPU,
stills and rotation), and derive each from its own axis.

On a 149/83/226 A cell that is 23.1M candidates per frame instead of 94.2M, 4.1x
fewer. Results are bit-identical, as they must be - the candidates removed are only
ones the |q| <= 1/d_min cut rejected anyway: over six rotation crystals every merged
observation count, high-shell CC1/2 and space group matches the cube exactly, 6/6
space groups correct.

It buys almost no time, and the earlier claim that the cube cost 22% of that
crystal's wall clock was wrong. Removing 4.1x of the candidates moves it 1m58s ->
1m57s, so the whole prediction sweep is ~1% of the run. The 22% that crystal costs
relative to a fixed max_hkl of 100 is genuine extra work at max_l = 227: real
reflections inside the resolution sphere along the long axis, predicted and
integrated either way. Per-axis limits do not reduce that and cannot.

The user-facing setting stays a single number: it exists to bound the work, not to
describe the crystal, and applies to all three indices when set.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

image_analysis/IndexAndRefine.cpp

@@ -22,29 +22,32 @@ namespace {
     // so this bounds the pathological case rather than the normal one.
     constexpr int OFFLINE_REFINE_ITERATIONS = 50;
 
-    // How far the predictor has to walk the lattice for THIS crystal. The predictor keeps only
+    // How far the predictor has to walk each index for THIS crystal. The predictor keeps only
     // reflections with |q| <= 1/d_min, and h = a.q for the real-space axis a, so |h| <= a/d_min exactly
-    // (Cauchy-Schwarz, equality when a lies along q); the same for k and l. Walking beyond
-    // max(a,b,c)/d_min therefore generates only candidates the resolution cut throws away, and stopping
-    // short of it silently truncates the outermost reflections of the longest axis.
-    //
-    // A fixed bound cannot be right for both: 100 covers a 150 A axis at 1.5 A and truncates the same
-    // axis at 1.0 A, while being ~4x more cube than a small cell at 2 A ever needs. One index of margin
-    // covers the rounding.
-    int MaxHKLForCell(const CrystalLattice &latt, float d_min_A) {
-        const UnitCell cell = latt.GetUnitCell();
-        const float longest_axis_A = std::max({cell.a, cell.b, cell.c});
-        return static_cast<int>(std::ceil(longest_axis_A / d_min_A)) + 1;
+    // (Cauchy-Schwarz, equality when a lies along q) - and independently |k| <= b/d_min, |l| <= c/d_min.
+    // Each index is therefore bounded by its OWN axis, which is why the limits are per-axis: a single
+    // cube would have to be sized for the longest axis and would walk the short ones far past anything
+    // the resolution cut can keep. One index of margin covers the rounding.
+    int MaxIndexForAxis(float axis_A, float d_min_A) {
+        return static_cast<int>(std::ceil(axis_A / d_min_A)) + 1;
     }
 
-    // An explicit setting is enforced as given; otherwise the cell decides. The online path carries a
-    // value (the broker bootstraps one and the API can change it), so live acquisition never has its
-    // per-frame cost decided by the crystal that happened to be mounted.
-    int PredictionMaxHKL(const DiffractionExperiment &experiment, const CrystalLattice &latt) {
+    // An explicit setting is enforced as given, on every index - it is one number, deliberately, because
+    // it exists to bound the work rather than to describe the crystal. Otherwise the cell decides. The
+    // online path always carries a value (the broker bootstraps one and the API can change it), so a
+    // live acquisition never has its per-frame cost decided by whichever crystal was mounted.
+    void ApplyPredictionRange(BraggPredictionSettings &settings, const DiffractionExperiment &experiment,
+                              const CrystalLattice &latt) {
         const auto &bragg = experiment.GetBraggIntegrationSettings();
-        if (const auto fixed = bragg.GetMaxHKL())
-            return *fixed;
-        return MaxHKLForCell(latt, bragg.GetDMinLimit_A());
+        if (const auto fixed = bragg.GetMaxHKL()) {
+            settings.max_h = settings.max_k = settings.max_l = *fixed;
+            return;
+        }
+        const UnitCell cell = latt.GetUnitCell();
+        const float d_min_A = bragg.GetDMinLimit_A();
+        settings.max_h = MaxIndexForAxis(cell.a, d_min_A);
+        settings.max_k = MaxIndexForAxis(cell.b, d_min_A);
+        settings.max_l = MaxIndexForAxis(cell.c, d_min_A);
     }
 }
 
@@ -462,10 +465,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
         .image_scale_cc = msg.image_scale_cc,
     };
 
-    const BraggPredictionSettings settings_prediction{
+    BraggPredictionSettings settings_prediction{
         .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
         .ewald_dist_cutoff = ewald_dist_cutoff,
-        .max_hkl = PredictionMaxHKL(experiment, latt),
         // Centering is a hypothesis to confirm, not assume: with no user-fixed space group, predict
         // in P so the centering-absent reflections are integrated and the space-group search can
         // confirm or disprove centering (and catch a missed superstructure). A user-fixed space
@@ -476,6 +478,7 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
         // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged.
         .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f,
     };
+    ApplyPredictionRange(settings_prediction, experiment, latt);
 
     // Predict, then integrate with the selected integrator (box-sum or profile-fit).
     auto pred_start_time = std::chrono::steady_clock::now();