RESULT-CHANGING (owner approved, ulp level): CPU profile fit in the GPU engine's arithmetic

Pass B of BraggIntegrationEngineCPU now does what BraggIntegrationEngineGPU's fit does: profiles in
float, the Gaussian from expf with the widths' reciprocals taken once and normalised by a reciprocal,
and the reweighted fit in float with one reciprocal of the variance for its three sums (the GPU's
__frcp_rn). Per-reflection geometry and the final sigma stay in double. The logic was already the
same on both; the CPU now differs from the GPU only in summation order.

Measured (CPU-only, Apple M6): Pass B loops 1.2x faster on P cores; cytochrome C 1.9%, myoglobin 1.5%
faster end to end, 4nwv integration steps 3%. Intensities move by ~1e-6 relative; merged statistics
in the fourth digit (error model b 1.3259e-3 -> 1.3258e-3, R_work 0.2440 -> 0.2439).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -537,22 +537,22 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     };
     // Normalised empirical profile (sum = 1), the average grid over the strong spots of a shell.
     // ProfileGaussian does not use it - it rebuilds a per-reflection Gaussian in Pass B.
+    // Stored in float, as the fit below reads it (and as the GPU engine keeps it).
     auto build_profile = [&](const std::vector<double> &grid) {
-        std::vector<double> P(GG, 0.0);
+        std::vector<float> P(GG, 0.0f);
         double sum = 0.0;
-        for (int k = 0; k < GG; ++k) {
-            P[k] = std::max(0.0, grid[k]);
-            sum += P[k];
-        }
+        for (int k = 0; k < GG; ++k)
+            sum += std::max(0.0, grid[k]);
         if (sum > 0.0)
-            for (double &p : P) p /= sum;
+            for (int k = 0; k < GG; ++k)
+                P[k] = static_cast<float>(std::max(0.0, grid[k]) / sum);
         return P;
     };
 
-    const std::vector<double> global_P = empirical && global_n > 0 ? build_profile(global_grid)
-                                                                   : std::vector<double>(GG, 0.0);
+    const std::vector<float> global_P = empirical && global_n > 0 ? build_profile(global_grid)
+                                                                  : std::vector<float>(GG, 0.0f);
     const Sigma2 global_sigma2 = widths(global_mom);
-    std::vector<std::vector<double>> shell_P(N_SHELL, global_P);
+    std::vector<std::vector<float>> shell_P(N_SHELL, global_P);
     std::vector<Sigma2> shell_sigma2(N_SHELL, global_sigma2);
     for (int s = 0; s < N_SHELL; ++s) {
         if (shell_n[s] < MIN_STRONG_PER_SHELL) continue;
@@ -562,11 +562,14 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
 
     // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). The
     //     reweighting is the Kabsch/Otwinowski iteration: Kabsch, Acta Cryst D66, 133-144 (2010);
-    //     Otwinowski & Minor, Methods Enzymol 276, 307-326 (1997). ---
-    std::vector<double> Pbuf;
+    //     Otwinowski & Minor, Methods Enzymol 276, 307-326 (1997).
+    //     In single precision, with the GPU engine's arithmetic (BraggIntegrationEngineGPU.cu, fit):
+    //     the profile from expf with the widths' reciprocals taken once, normalised by a reciprocal,
+    //     and one reciprocal of the variance for the three weights. ---
+    std::vector<float> Pbuf;
     // The pixels the fit reads, gathered once while p_valid walks the grid: the profile value and the
     // background-subtracted count, in the grid's order.
-    struct FitPixel { double Pp, signal; };
+    struct FitPixel { float Pp, signal; };
     std::vector<FitPixel> fit_pixels;
     for (size_t i = 0; i < npredicted; ++i) {
         const auto &rh = rough[i];
@@ -584,7 +587,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         const int sh = rh.shell < 0 ? 0 : rh.shell;
 
         int Rf = R;
-        const std::vector<double> *Pvec = &shell_P[sh];
+        const std::vector<float> *Pvec = &shell_P[sh];
         const BraggStencil st = MakeBraggStencil(predicted[i].predicted_x, predicted[i].predicted_y, stencil);
         if (!empirical) {
             const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
@@ -625,17 +628,21 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                 prefetch_row(img.address(static_cast<size_t>(y) * W + gx0),
                              img.address(static_cast<size_t>(y) * W + gx1));
             const int Gf = 2 * Rf + 1;
-            Pbuf.assign(static_cast<size_t>(Gf) * Gf, 0.0);
-            double gs = 0.0;
+            Pbuf.assign(static_cast<size_t>(Gf) * Gf, 0.0f);
+            const float fxf = static_cast<float>(fx), fyf = static_cast<float>(fy);
+            const float uxf = static_cast<float>(ux), uyf = static_cast<float>(uy);
+            const float inv_2s2r = 0.5f / static_cast<float>(s2r), inv_2s2t = 0.5f / static_cast<float>(s2t);
+            float gs = 0.0f;
             for (int dy = -Rf; dy <= Rf; ++dy)
                 for (int dx = -Rf; dx <= Rf; ++dx) {
-                    const double ex = dx - fx, ey = dy - fy;
-                    const double rad = ex * ux + ey * uy, tn = -ex * uy + ey * ux;
-                    const double g = std::exp(-rad * rad / (2.0 * s2r) - tn * tn / (2.0 * s2t));
+                    const float ex = dx - fxf, ey = dy - fyf;
+                    const float rad = ex * uxf + ey * uyf, tn = -ex * uyf + ey * uxf;
+                    const float g = std::exp(-rad * rad * inv_2s2r - tn * tn * inv_2s2t);
                     Pbuf[(dy + Rf) * Gf + (dx + Rf)] = g;
                     gs += g;
                 }
-            for (double &p : Pbuf) p /= gs;
+            const float inv_gs = 1.0f / gs;
+            for (float &p : Pbuf) p *= inv_gs;
             Pvec = &Pbuf;
         }
 
@@ -648,13 +655,13 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         //     pixels, what Reject cuts on. m_read / m_all is the fraction of the r1 disk the
         //     summation seed the runaway guard compares against actually summed; with nothing missing
         //     and nothing excluded it is 1 and the guard is untouched. ---
-        double p_grid = 0.0, p_valid = 0.0, p_own = 0.0, m_all = 0.0, m_read = 0.0;
-        double p_peak = 0.0, p_lost_peak = 0.0;
+        float p_grid = 0.0f, p_valid = 0.0f, p_own = 0.0f, m_all = 0.0f, m_read = 0.0f;
+        float p_peak = 0.0f, p_lost_peak = 0.0f;
         fit_pixels.clear();
         for (int dy = -Rf; dy <= Rf; ++dy)
             for (int dx = -Rf; dx <= Rf; ++dx) {
-                const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)];
-                if (Pp <= 0.0) continue;
+                const float Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)];
+                if (Pp <= 0.0f) continue;
                 const int x = rh.cx + dx, y = rh.cy + dy;
                 if (x < 0 || y < 0 || x >= W || y >= H) continue;
                 // The region the summation seed was taken over (pass A), so the guard below
@@ -677,7 +684,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                 if (own) p_own += Pp;
                 if (in_disk && (own || !exclude)) m_read += Pp;
                 if (own || !exclude)
-                    fit_pixels.push_back({Pp, static_cast<double>(img[y * W + x]) - rh.bkg});
+                    fit_pixels.push_back({Pp, static_cast<float>(img[y * W + x]) - static_cast<float>(rh.bkg)});
             }
         if (p_valid < overlap_min_peak * p_grid) continue;
         // A hole in the profile's PEAK is a different defect from a hole in its wings, and the mass
@@ -687,21 +694,23 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         if (p_lost_peak > MINPK_MAX_MISSING_PEAK * p_peak) continue;
         if (overlap == OverlapMode::Reject && p_own < overlap_min_peak) continue;
 
-        const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR);
-        double I = rh.I, den = 0.0, wsum = 0.0;
+        const float B = static_cast<float>(std::max(rh.bkg, PIXEL_VARIANCE_FLOOR));
+        const float v_min = static_cast<float>(WEIGHT_VARIANCE_MIN_FRACTION) * B;
+        float If = static_cast<float>(rh.I), den = 0.0f, wsum = 0.0f;
         for (int iter = 0; iter < 4; ++iter) {
-            double num = 0.0;
-            den = 0.0;
-            wsum = 0.0;
+            float num = 0.0f;
+            den = 0.0f;
+            wsum = 0.0f;
             for (const auto &[Pp, signal] : fit_pixels) {
-                const double v = std::max(B + I * Pp, WEIGHT_VARIANCE_MIN_FRACTION * B);
-                num += Pp * signal / v;
-                den += Pp * Pp / v;
-                wsum += Pp / v;
+                const float iv = 1.0f / std::max(B + If * Pp, v_min);
+                num += Pp * signal * iv;
+                den += Pp * Pp * iv;
+                wsum += Pp * iv;
             }
-            if (den > 0.0) I = num / den; else break;
+            if (den > 0.0f) If = num / den; else break;
         }
-        if (!(den > 0.0)) continue;
+        if (!(den > 0.0f)) continue;
+        double I = If;
 
         // Guard against profile-fit runaways: on a weak / near-zero reflection the reweighted Kabsch
         // iteration has no real peak to lock onto and can manufacture intensity the box sum never sees.