RESULT-CHANGING (owner approved, ulp level): Metal profile fit, CPU Gaussian from a shared exp

Pass B of BraggIntegrationEngineCPU is split where a device can take over: per reflection, what the
fit needs is worked out in double (BraggFitSetup), the profile, its readable mass and the reweighted
fit run in single precision (FitProfile -> BraggFitOut), and the error, the summation guard and the
result are finished in double. BraggIntegrationEngineMetal runs the middle as a Metal kernel, one
thread per reflection, summing in the CPU's order (FitOnDevice); the CUDA engine is untouched.

The Gaussian's exp is now BraggExp (BraggExp.h), on the CPU and on the GPU: no GPU exp reproduces the
host's expf, and BraggExp uses only +, -, * and exact power-of-two scaling, so both compile it to the
same bits. Checked on the Apple M6 over every float in [-110, 89] (2.2e9 arguments), with subnormal
products and sums matching too. It is closer to correctly rounded than macOS expf (max 0.53 ulp;
0.2% of arguments below -1 not correctly rounded vs 0.62%). This is the one result change: per
observation |dI|/sigma p99.9 <= 1.7e-6, max 1e-4, 10-30x below what ba3f75475 moved. Merged
statistics (resolution, completeness, CC1/2, R_meas, I/sigma, ISa, space group) are identical on
lysozyme, cytochrome C, myoglobin, a small-molecule set and 4nwv; only tNCS fields that flip on any
last-digit change move (cytochrome C's pseudo-translation flag, whose modulation ratio sits at the
2.5 threshold, turns off).

The Metal engine stays bit-identical to the CPU engine: BraggIntegrationEngineMetalTest passes
unchanged, and catches planted kernel bugs (another exp, 3 iterations, a dropped lost-peak test, <=
for the disk, inverted exclude logic).

Measured wall time, Apple M6 (CPU pass B -> GPU pass B): lysozyme 24.7 -> 18.6 s, myoglobin 39.8 ->
35.7 s, cytochrome C 59.3 -> 49.0 s, 4nwv 150.2 -> 131.3 s; the small-molecule set, sparse frames,
92.5 -> 95.6 s (noise of +-3 s). Metal's fast::exp in the kernel instead measured the same on every
set (18.6, 35.8, 49.1, 131.2 s), so it would buy nothing for the bitwise check it would cost. On the
CPU alone BraggExp costs nothing measurable (cytochrome C --no-gpu 72.8 -> 73.1 s).

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

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -429,6 +429,107 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection
     a.rough[i] = out;
 }
 
+bool BraggIntegrationEngineCPU::FitKept(const BraggFitOut &o) const {
+    if (o.p_valid < overlap_min_peak * o.p_grid) return false;
+    // A hole in the profile's PEAK is a different defect from a hole in its wings, and the mass
+    // fraction above cannot tell them apart - the peak of a broad spot is a few percent of the
+    // mass, so MINPK passes a reflection that has lost the one part of the profile its amplitude
+    // is determined by. See the header.
+    if (o.p_lost_peak > MINPK_MAX_MISSING_PEAK * o.p_peak) return false;
+    if (overlap == OverlapMode::Reject && o.p_own < overlap_min_peak) return false;
+    return true;
+}
+
+template <class Sampler, class OwnerMap>
+BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const OwnerMap &owner_map,
+                                                  const BraggFitSetup &s, size_t i,
+                                                  const std::vector<float> &shell_profiles,
+                                                  std::vector<float> &cells,
+                                                  std::vector<FitPixel> &fit_pixels) const {
+    const int W = static_cast<int>(xpixel), H = static_cast<int>(ypixel);
+    const bool exclude = overlap == OverlapMode::Exclude;
+    const bool use_owner = overlap != OverlapMode::Off;
+    const int Rf = s.Rf, Gf = 2 * Rf + 1;
+    BraggFitOut o{};
+
+    const float *P;
+    if (s.shell >= 0)
+        P = shell_profiles.data() + static_cast<size_t>(s.shell) * GG;
+    else {
+        // The grid's image rows, asked for now (as in pass A): building the profile below takes long
+        // enough for them to arrive before the fit reads them.
+        const int gx0 = std::max(0, s.cx - Rf), gx1 = std::min(W - 1, s.cx + Rf);
+        for (int y = std::max(0, s.cy - Rf); y <= std::min(H - 1, s.cy + Rf); ++y)
+            prefetch_row(img.address(static_cast<size_t>(y) * W + gx0),
+                         img.address(static_cast<size_t>(y) * W + gx1));
+        cells.resize(static_cast<size_t>(Gf) * Gf);
+        float gs = 0.0f;
+        for (int dy = -Rf; dy <= Rf; ++dy)
+            for (int dx = -Rf; dx <= Rf; ++dx) {
+                const float g = BraggGaussian(dx, dy, s.fx, s.fy, s.ux, s.uy, s.inv_2s2r, s.inv_2s2t);
+                cells[(dy + Rf) * Gf + (dx + Rf)] = g;
+                gs += g;
+            }
+        const float inv_gs = 1.0f / gs;
+        for (float &p : cells) p *= inv_gs;
+        P = cells.data();
+    }
+
+    // --- How much of the expected profile the fit can actually see. p_valid is the readable
+    //     fraction of it - XDS's MINPK, dials' valid_foreground_threshold - measured against the
+    //     mass that falls on the detector at all, so a reflection is judged on the pixels that
+    //     exist and not on where the sensor ends. p_own is the same quantity over neighbour-owned
+    //     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. ---
+    const BraggStencil st = MakeBraggStencil(s.px, s.py, stencil);
+    fit_pixels.clear();
+    for (int dy = -Rf; dy <= Rf; ++dy)
+        for (int dx = -Rf; dx <= Rf; ++dx) {
+            const float Pp = P[(dy + Rf) * Gf + (dx + Rf)];
+            if (Pp <= 0.0f) continue;
+            const int x = s.cx + dx, y = s.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 compares the two
+            // over the same pixels.
+            const bool in_disk = s.footprint
+                    ? BraggStencilInSignal(st, x - s.px, y - s.py,
+                                           (x - s.px) * (x - s.px) + (y - s.py) * (y - s.py), r1_sq)
+                    : dx * dx + dy * dy < r1_sq;
+            o.p_grid += Pp;
+            o.p_peak = std::max(o.p_peak, Pp);
+            if (in_disk) o.m_all += Pp;
+            const int32_t v = img[y * W + x];
+            if (!valid(v)) {
+                o.p_lost_peak = std::max(o.p_lost_peak, Pp);
+                continue;
+            }
+            o.p_valid += Pp;
+            const bool own = !use_owner || BraggOwnedBy(owner_map.Get(x, y), static_cast<int>(i));
+            if (own) o.p_own += Pp;
+            if (in_disk && (own || !exclude)) o.m_read += Pp;
+            if (own || !exclude)
+                fit_pixels.push_back({Pp, static_cast<float>(v) - s.bkg});
+        }
+    if (!FitKept(o)) return o;
+
+    float If = s.I0;
+    for (int iter = 0; iter < 4; ++iter) {
+        float num = 0.0f;
+        o.den = 0.0f;
+        o.wsum = 0.0f;
+        for (const auto &[Pp, signal] : fit_pixels) {
+            const float iv = 1.0f / std::max(s.B + If * Pp, s.v_min);
+            num += Pp * signal * iv;
+            o.den += Pp * Pp * iv;
+            o.wsum += Pp * iv;
+        }
+        if (o.den > 0.0f) If = num / o.den; else break;
+    }
+    o.I = If;
+    return o;
+}
+
 template <class Sampler, class OwnerMap>
 std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &img, const OwnerMap &owner_map,
                                                                  const std::vector<Reflection> &predicted,
@@ -578,16 +679,21 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
     // --- 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).
-    //     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 { float Pp, signal; };
-    std::vector<FitPixel> fit_pixels;
+    //     Per reflection in three steps: what the fit needs, worked out in double (BraggFitSetup); the
+    //     fit in single precision (FitProfile, or the Metal engine's kernel), with the GPU engine's
+    //     arithmetic (BraggIntegrationEngineGPU.cu, fit) - the profile from exp with the widths'
+    //     reciprocals taken once, normalised by a reciprocal, and one reciprocal of the variance for the
+    //     three weights; and its error, the summation guard and the result, in double. ---
+    std::vector<float> shell_profiles;   // the learned profiles, shell by shell (ProfileEmpirical)
+    if (empirical)
+        for (int s = 0; s < N_SHELL; ++s)
+            shell_profiles.insert(shell_profiles.end(), shell_P[s].begin(), shell_P[s].end());
+    std::vector<BraggFitSetup> setups(npredicted);
+    uint32_t cells = 0;
     for (size_t i = 0; i < npredicted; ++i) {
         const auto &rh = rough[i];
+        BraggFitSetup &fs = setups[i];
+        fs = {};
         if (!rh.ok) continue;
         // An overloaded reflection is not measured by any fit - its brightest pixels are what is
         // missing. It is passed on with its box sum and the flag, for the merge to drop its rocking
@@ -600,10 +706,15 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
             continue;
         }
         const int sh = rh.shell < 0 ? 0 : rh.shell;
-
-        int Rf = R;
-        const std::vector<float> *Pvec = &shell_P[sh];
         const BraggStencil st = MakeBraggStencil(predicted[i].predicted_x, predicted[i].predicted_y, stencil);
+        fs.fit = 1;
+        fs.px = predicted[i].predicted_x;
+        fs.py = predicted[i].predicted_y;
+        fs.cx = rh.cx;
+        fs.cy = rh.cy;
+        fs.Rf = R;
+        fs.shell = sh;
+        fs.footprint = st.fp_s2t > 0.0f;
         if (!empirical) {
             const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
             const double Rpx = std::hypot(rx, ry);
@@ -632,100 +743,41 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
                 s2r = std::max<double>(s2r, st.fp_s2r);
                 elong = true;
             }
-            // Build the Gaussian per reflection, centred on the sub-pixel predicted position and (when
+            // The Gaussian is built per reflection, centred on the sub-pixel predicted position and (when
             // needed) elongated, on a grid grown to hold the wider of its two axes.
-            const double fx = predicted[i].predicted_x - rh.cx, fy = predicted[i].predicted_y - rh.cy;
-            Rf = elong ? std::min(3 * R, static_cast<int>(std::ceil(r2 + 2.0 * std::sqrt(std::max(s2r, s2t))))) : R;
-            // The grid's image rows, asked for now (as in pass A): building the profile below takes
-            // long enough for them to arrive before the fit reads them.
-            const int gx0 = std::max(0, rh.cx - Rf), gx1 = std::min(W - 1, rh.cx + Rf);
-            for (int y = std::max(0, rh.cy - Rf); y <= std::min(H - 1, rh.cy + Rf); ++y)
-                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.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 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;
-                }
-            const float inv_gs = 1.0f / gs;
-            for (float &p : Pbuf) p *= inv_gs;
-            Pvec = &Pbuf;
+            fs.Rf = elong ? std::min(3 * R, static_cast<int>(std::ceil(r2 + 2.0 * std::sqrt(std::max(s2r, s2t))))) : R;
+            fs.shell = -1;
+            fs.cells = cells;
+            cells += (2 * fs.Rf + 1) * (2 * fs.Rf + 1);
+            fs.fx = static_cast<float>(predicted[i].predicted_x - rh.cx);
+            fs.fy = static_cast<float>(predicted[i].predicted_y - rh.cy);
+            fs.ux = static_cast<float>(ux);
+            fs.uy = static_cast<float>(uy);
+            fs.inv_2s2r = 0.5f / static_cast<float>(s2r);
+            fs.inv_2s2t = 0.5f / static_cast<float>(s2t);
         }
+        fs.bkg = static_cast<float>(rh.bkg);
+        fs.B = static_cast<float>(std::max(rh.bkg, PIXEL_VARIANCE_FLOOR));
+        fs.v_min = static_cast<float>(WEIGHT_VARIANCE_MIN_FRACTION) * fs.B;
+        fs.I0 = static_cast<float>(rh.I);
+    }
 
-        const int Gf = 2 * Rf + 1;
+    std::vector<BraggFitOut> fits(npredicted);
+    if (!FitOnDevice(setups, shell_profiles, fits)) {
+        std::vector<float> cells_buf;
+        std::vector<FitPixel> fit_pixels;
+        for (size_t i = 0; i < npredicted; ++i)
+            if (setups[i].fit)
+                fits[i] = FitProfile(img, owner_map, setups[i], i, shell_profiles, cells_buf, fit_pixels);
+    }
 
-        // --- How much of the expected profile the fit can actually see. p_valid is the readable
-        //     fraction of it - XDS's MINPK, dials' valid_foreground_threshold - measured against the
-        //     mass that falls on the detector at all, so a reflection is judged on the pixels that
-        //     exist and not on where the sensor ends. p_own is the same quantity over neighbour-owned
-        //     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. ---
-        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 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
-                // compares the two over the same pixels.
-                const bool in_disk = st.fp_s2t > 0.0f
-                        ? BraggStencilInSignal(st, x - predicted[i].predicted_x, y - predicted[i].predicted_y,
-                                               (x - predicted[i].predicted_x) * (x - predicted[i].predicted_x)
-                                               + (y - predicted[i].predicted_y) * (y - predicted[i].predicted_y),
-                                               r1_sq)
-                        : dx * dx + dy * dy < r1_sq;
-                p_grid += Pp;
-                p_peak = std::max(p_peak, Pp);
-                if (in_disk) m_all += Pp;
-                if (!valid(img[y * W + x])) {
-                    p_lost_peak = std::max(p_lost_peak, Pp);
-                    continue;
-                }
-                p_valid += Pp;
-                const bool own = clean(x, y, i);
-                if (own) p_own += Pp;
-                if (in_disk && (own || !exclude)) m_read += Pp;
-                if (own || !exclude)
-                    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
-        // fraction above cannot tell them apart - the peak of a broad spot is a few percent of the
-        // mass, so MINPK passes a reflection that has lost the one part of the profile its amplitude
-        // is determined by. See the header.
-        if (p_lost_peak > MINPK_MAX_MISSING_PEAK * p_peak) continue;
-        if (overlap == OverlapMode::Reject && p_own < overlap_min_peak) continue;
-
-        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) {
-            float num = 0.0f;
-            den = 0.0f;
-            wsum = 0.0f;
-            for (const auto &[Pp, signal] : fit_pixels) {
-                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.0f) If = num / den; else break;
-        }
-        if (!(den > 0.0f)) continue;
-        double I = If;
+    for (size_t i = 0; i < npredicted; ++i) {
+        if (!setups[i].fit) continue;
+        const auto &rh = rough[i];
+        const BraggFitOut &fo = fits[i];
+        if (!FitKept(fo) || !(fo.den > 0.0f)) continue;
+        const float den = fo.den, wsum = fo.wsum;
+        double I = fo.I;
 
         // 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.
@@ -746,7 +798,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
         // The seed is a sum over the disk the box sum actually read, so when pixels are missing from
         // both - excluded to a neighbour, or unreadable - the fit's full-profile intensity has to be
         // scaled down to that same disk before the two are comparable. Nothing dropped gives 1.
-        const double guard_scale = m_all > 0.0 ? m_read / m_all : 1.0;
+        const double guard_scale = fo.m_all > 0.0f ? fo.m_read / fo.m_all : 1.0;
         if (std::abs(I * guard_scale - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) {
             I = rh.I;
             sigma = rh.sigma;