From ba3f7547520365622bddc91d2308594a7fde886f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 14 Oct 2026 19:42:31 +0200 Subject: [PATCH] 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 --- .../BraggIntegrationEngineCPU.cpp | 79 +++++++++++-------- 1 file changed, 44 insertions(+), 35 deletions(-) diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index b62ab2aa9..4b7f785a9 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -537,22 +537,22 @@ std::vector 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 &grid) { - std::vector P(GG, 0.0); + std::vector 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(std::max(0.0, grid[k]) / sum); return P; }; - const std::vector global_P = empirical && global_n > 0 ? build_profile(global_grid) - : std::vector(GG, 0.0); + const std::vector global_P = empirical && global_n > 0 ? build_profile(global_grid) + : std::vector(GG, 0.0f); const Sigma2 global_sigma2 = widths(global_mom); - std::vector> shell_P(N_SHELL, global_P); + std::vector> shell_P(N_SHELL, global_P); std::vector 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 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 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 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 fit_pixels; for (size_t i = 0; i < npredicted; ++i) { const auto &rh = rough[i]; @@ -584,7 +587,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const int sh = rh.shell < 0 ? 0 : rh.shell; int Rf = R; - const std::vector *Pvec = &shell_P[sh]; + const std::vector *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 BraggIntegrationEngineCPU::RunImpl(const Sampler &img, prefetch_row(img.address(static_cast(y) * W + gx0), img.address(static_cast(y) * W + gx1)); const int Gf = 2 * Rf + 1; - Pbuf.assign(static_cast(Gf) * Gf, 0.0); - double gs = 0.0; + Pbuf.assign(static_cast(Gf) * Gf, 0.0f); + const float fxf = static_cast(fx), fyf = static_cast(fy); + const float uxf = static_cast(ux), uyf = static_cast(uy); + const float inv_2s2r = 0.5f / static_cast(s2r), inv_2s2t = 0.5f / static_cast(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 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 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(img[y * W + x]) - rh.bkg}); + fit_pixels.push_back({Pp, static_cast(img[y * W + x]) - static_cast(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 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(std::max(rh.bkg, PIXEL_VARIANCE_FLOOR)); + const float v_min = static_cast(WEIGHT_VARIANCE_MIN_FRACTION) * B; + float If = static_cast(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.