Accelerated Integration

CMakeLists.txt

@@ -8,6 +8,28 @@ include(ResolveGitVersion)
 
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CUDA_STANDARD 20)
+
+# Set CUDA compute capability target
+# Override with just the number: cmake -DCUDA_ARCH=86 ..
+# Set default to 7.5 (Turing)
+set(CUDA_ARCH "75 - Turing (RTX 20xx, GTX 16xx)" CACHE STRING "CUDA compute capability target architecture")
+set_property(CACHE CUDA_ARCH PROPERTY STRINGS 
+    "52 - Maxwell (GTX 9xx, Titan X)"
+    "60 - Pascal (GTX 10xx, P100)"
+    "61 - Pascal (GTX 1050/1030)"
+    "70 - Volta (V100, Titan V100)"
+    "75 - Turing (RTX 20xx, GTX 16xx)"
+    "80 - Ampere (A100)"
+    "86 - Ampere (RTX 30xx)"
+    "89 - Ada Lovelace (RTX 40xx)"
+    "90 - Hopper (H100)" 
+)
+
+# Extract just the number from the selection (first two characters)
+string(SUBSTRING "${CUDA_ARCH}" 0 2 CUDA_ARCH_NUMBER)
+set(CMAKE_CUDA_ARCHITECTURES ${CUDA_ARCH_NUMBER})
+message(STATUS "CUDA architecture set to: sm_${CMAKE_CUDA_ARCHITECTURES}")
+
 project(fast-feedback-service LANGUAGES CXX VERSION ${FFS_VERSION_CMAKE})
 
 set(CMAKE_EXPORT_COMPILE_COMMANDS yes)
@@ -111,6 +133,36 @@ find_package(CUDAToolkit)
 if (CUDAToolkit_FOUND)
     message(STATUS "CUDA found: Building CUDA components.")
     set(FFS_ENABLE_CUDA ON)
+
+    # Try to detect GPU compute capability
+    execute_process(
+        COMMAND nvidia-smi --query-gpu=compute_cap --format=csv,noheader
+        OUTPUT_VARIABLE DETECTED_GPU_CAPS
+        OUTPUT_STRIP_TRAILING_WHITESPACE
+        ERROR_QUIET
+    )
+
+    if(DETECTED_GPU_CAPS)
+        # nvidia-smi returns "7.5" format, convert to "75"
+        string(REPLACE "." "" DETECTED_GPU_ARCH "${DETECTED_GPU_CAPS}")
+        # Get first GPU if multiple
+        string(REGEX MATCH "[0-9]+" DETECTED_GPU_ARCH "${DETECTED_GPU_ARCH}")
+        message(STATUS "Detected GPU compute capability: sm_${DETECTED_GPU_ARCH}")
+
+        # Warn if configured architecture doesn't match detected GPU
+        if(NOT CUDA_ARCH_NUMBER STREQUAL DETECTED_GPU_ARCH)
+            message(WARNING 
+                "Target CUDA architecture (sm_${CUDA_ARCH_NUMBER}) differs from detected GPU (sm_${DETECTED_GPU_ARCH}).\n"
+                "   The compiled code may not run optimally or at all on this system.\n"
+                "   To build for the detected GPU: cmake -DCUDA_ARCH=${DETECTED_GPU_ARCH} ..\n"
+                "   To build for a different target: ensure the target system has sm_${CUDA_ARCH_NUMBER} or newer."
+            )
+        else()
+            message(STATUS "✓ Configured architecture matches detected GPU")
+        endif()
+    else()
+        message(STATUS "Could not detect GPU compute capability (nvidia-smi not available or no GPU present)")
+    endif()
 else()
     message(STATUS "CUDA not found: Skipping CUDA components.")
     set(FFS_ENABLE_CUDA OFF)

include/math/device_precision.cuh

@@ -34,4 +34,15 @@ using scalar_t = float;
 #define CUDA_ABS fabsf
 #define CUDA_SQRT sqrtf
 #define CUDA_EXP expf
-#endif
+#endif
+
+// Accumulator type for summation/reduction operations.
+// Integer accumulation is used because:
+// 1. atomicAdd on uint32_t is a single native hardware instruction (fastest atomic)
+// 2. Maximum possible sum is well within uint32_t range (~20M << 4.3G)
+// 3. Final reduction to double is performed on the host
+//
+// NOTE: If pixel_t can be negative (pedestal-subtracted data), change to int32_t.
+// NOTE: If pixel values are non-integer (e.g. gain-corrected floats), revert to:
+//     using accumulator_t = scalar_t;
+using accumulator_t = uint32_t;

include/math/math_utils.cuh

@@ -13,6 +13,7 @@
  * constexpr scalar_t RAD_TO_DEG = 180.0 / cuda::std::numbers::pi_v<scalar_t>;
  */
 #include <cmath>
+#include <concepts>
 
 #ifdef __CUDACC__
 #include "math/device_precision.cuh"
@@ -43,4 +44,24 @@ DEVICE_HOST scalar_type degrees_to_radians(scalar_type degrees) {
     return degrees * DEG_TO_RAD;
 }
 
+/**
+ * @brief Integer ceiling division
+ * 
+ * Computes ceil(n/d) using integer arithmetic only, avoiding
+ * floating-point conversion. Equivalent to (n + d - 1) / d, which
+ * rounds up to the nearest integer quotient.
+ * 
+ * @tparam T Integer type
+ * @param n Numerator
+ * @param d Denominator (must be > 0)
+ * @return Ceiling of n/d
+ * 
+ * @example ceil_div(15, 4) returns 4, whereas 15/4 returns 3
+ */
+template <typename T>
+DEVICE_HOST constexpr T ceil_div(T n, T d) {
+    static_assert(std::is_integral_v<T>, "ceil_div requires an integer type");
+    return (n + d - 1) / d;
+}
+
 #undef DEVICE_HOST

integrator/CMakeLists.txt

@@ -6,6 +6,9 @@ find_package(Bitshuffle REQUIRED)
 
 add_executable(integrator
     integrator.cc
+    integrator.cu
+    kabsch.cu
+    extent.cc
     ../spotfinder/shmread.cc
     ../spotfinder/cbfread.cc
 )
@@ -39,4 +42,7 @@ target_link_libraries(integrator
 # Allow CUDA to use multi-file device variables
 set_property(TARGET integrator PROPERTY CUDA_SEPARABLE_COMPILATION ON)
 target_compile_options(integrator PRIVATE "$<$<AND:$<CONFIG:Debug>,$<COMPILE_LANGUAGE:CUDA>>:-G>")
-target_compile_options(integrator PRIVATE "$<$<AND:$<COMPILE_LANGUAGE:CUDA>,$<OR:$<CONFIG:Debug>,$<CONFIG:RelWithDebInfo>>>:--generate-line-info>")
+target_compile_options(integrator PRIVATE "$<$<AND:$<COMPILE_LANGUAGE:CUDA>,$<OR:$<CONFIG:Debug>,$<CONFIG:RelWithDebInfo>>>:--generate-line-info>")
+
+# Add tests subdirectory
+add_subdirectory(tests)

integrator/extent.cc

@@ -0,0 +1,191 @@
+/**
+ * @file extent.cc
+ * @brief Extent and bounding box algorithms for baseline CPU implementation
+ */
+
+#include "extent.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+
+#include "ffs_logger.hpp"
+
+std::vector<BoundingBoxExtents> compute_kabsch_bounding_boxes(
+  const Eigen::Vector3d &s0,
+  const Eigen::Vector3d &rot_axis,
+  const mdspan_2d_double &s1_vectors,
+  const mdspan_2d_double &phi_positions,
+  const size_t num_reflections,
+  const double sigma_b,
+  const double sigma_m,
+  const Panel &panel,
+  const Scan &scan,
+  const MonochromaticBeam &beam,
+  const double n_sigma,
+  const double sigma_b_multiplier) {
+    std::vector<BoundingBoxExtents> extents;
+    extents.reserve(num_reflections);
+
+    /*
+    * Tolerance for detecting when a reflection is nearly parallel to
+    * the rotation axis. When ζ = m₂ · e₁ approaches zero, it indicates
+    * the reflection's scattering plane is nearly parallel to the
+    * goniometer rotation axis, making the φ-to-image conversion
+    * numerically unstable. This threshold (1e-10) is chosen based on
+    * geometric considerations rather than pure floating-point precision
+    * - it represents a practical limit for "nearly parallel" geometry
+    * where the standard bounding box calculation should be bypassed in
+    * favor of spanning the entire image range.
+    */
+    static constexpr double ZETA_TOLERANCE = 1e-10;
+
+    // Calculate the angular divergence parameters:
+    // Δb = nσ × σb × m (beam divergence extent)
+    // Δm = nσ × σm (mosaicity extent)
+    double delta_b = n_sigma * sigma_b * sigma_b_multiplier;
+    double delta_m = n_sigma * sigma_m;
+
+    // Extract experimental parameters needed for coordinate transformations
+    const auto [osc_start, osc_width] = scan.get_oscillation();
+    int image_range_start = scan.get_image_range()[0];
+    int image_range_end = scan.get_image_range()[1];
+    double wl = beam.get_wavelength();
+    Matrix3d d_matrix_inv = panel.get_d_matrix().inverse();
+
+    // Process each reflection individually
+    for (size_t i = 0; i < num_reflections; ++i) {
+        // Extract reflection centroid data
+        Eigen::Vector3d s1_c(
+          s1_vectors(i, 0), s1_vectors(i, 1), s1_vectors(i, 2));  // s₁ᶜ from s1_vectors
+        double phi_c = (phi_positions(i, 2));  // φᶜ from xyzcal.mm column
+
+        // Construct the Kabsch coordinate system for this reflection
+        // e1 = s₁ᶜ × s₀ / |s₁ᶜ × s₀| (perpendicular to scattering plane)
+        Eigen::Vector3d e1 = s1_c.cross(s0).normalized();
+        // e2 = s₁ᶜ × e₁ / |s₁ᶜ × e₁| (within scattering plane, orthogonal to e1)
+        Eigen::Vector3d e2 = s1_c.cross(e1).normalized();
+
+        double s1_len = s1_c.norm();
+
+        // Calculate s′ vectors at the four corners of the integration region
+        // These correspond to the extremes: (±Δb, ±Δb) in Kabsch coordinates
+        std::vector<Eigen::Vector3d> s_prime_vectors;
+        static constexpr std::array<std::pair<int, int>, 4> corner_signs = {
+          {{1, 1}, {1, -1}, {-1, 1}, {-1, -1}}};
+
+        for (auto [e1_sign, e2_sign] : corner_signs) {
+            // Project Δb divergences onto Kabsch basis vectors
+            // p represents the displacement in reciprocal space
+            Eigen::Vector3d p =
+              (e1_sign * delta_b * e1 / s1_len) + (e2_sign * delta_b * e2 / s1_len);
+
+            // Debug output for the Ewald sphere calculation
+            double p_magnitude = p.norm();
+            logger.trace(
+              "Reflection {}, corner ({},{}): p.norm()={:.6f}, s1_len={:.6f}, "
+              "delta_b={:.6f}",
+              i,
+              e1_sign,
+              e2_sign,
+              p_magnitude,
+              s1_len,
+              delta_b);
+
+            // Ensure the resulting s′ vector lies on the Ewald sphere
+            // This involves solving: |s′|² = |s₁ᶜ|² for the correct magnitude
+            double b = s1_len * s1_len - p.dot(p);
+            if (b < 0) {
+                logger.error(
+                  "Negative b value: {:.6f} for reflection {} (p.dot(p)={:.6f}, "
+                  "s1_len²={:.6f})",
+                  b,
+                  i,
+                  p.dot(p),
+                  s1_len * s1_len);
+                logger.error(
+                  "This means the displacement vector is too large for the Ewald "
+                  "sphere");
+                // Skip this corner or use a fallback approach
+                continue;
+            }
+            double d = -(p.dot(s1_c) / s1_len) + std::sqrt(b);
+
+            logger.trace("Reflection {}: b={:.6f}, d={:.6f}", i, b, d);
+
+            // Construct the s′ vector: s′ = (d × ŝ₁ᶜ) + p
+            Eigen::Vector3d s_prime = (d * s1_c / s1_len) + p;
+            s_prime_vectors.push_back(s_prime);
+        }
+
+        // Transform s′ vectors back to detector coordinates using Panel's get_ray_intersection
+        std::vector<std::pair<double, double>> detector_coords;
+        for (const auto &s_prime : s_prime_vectors) {
+            // Direct conversion from s′ vector to detector coordinates
+            // get_ray_intersection returns coordinates in mm
+            auto xy_mm_opt = panel.get_ray_intersection(s_prime);
+            if (!xy_mm_opt) {
+                continue;  // Skip this corner if no intersection
+            }
+            std::array<double, 2> xy_mm = *xy_mm_opt;
+
+            // Convert from mm to pixels using the new mm_to_px function
+            std::array<double, 2> xy_pixels = panel.mm_to_px(xy_mm[0], xy_mm[1]);
+
+            detector_coords.push_back({xy_pixels[0], xy_pixels[1]});
+        }
+
+        // Determine the bounding box in detector coordinates
+        // Find minimum and maximum coordinates from the four corners
+        auto [min_x_it, max_x_it] = std::minmax_element(
+          detector_coords.begin(),
+          detector_coords.end(),
+          [](const auto &a, const auto &b) { return a.first < b.first; });
+        auto [min_y_it, max_y_it] = std::minmax_element(
+          detector_coords.begin(),
+          detector_coords.end(),
+          [](const auto &a, const auto &b) { return a.second < b.second; });
+
+        BoundingBoxExtents bbox;
+        // Use floor/ceil as specified in the paper: xmin = floor(min([x1,x2,x3,x4]))
+        bbox.x_min = static_cast<int>(std::floor(min_x_it->first));
+        bbox.x_max = static_cast<int>(std::ceil(max_x_it->first));
+        bbox.y_min = static_cast<int>(std::floor(min_y_it->second));
+        bbox.y_max = static_cast<int>(std::ceil(max_y_it->second));
+
+        // Calculate the image range (z-direction) using mosaicity parameter Δm
+        // The extent in φ depends on the geometry factor ζ = m₂ · e₁
+        double zeta = rot_axis.dot(e1);
+        if (std::abs(zeta) > ZETA_TOLERANCE) {  // Avoid division by zero
+            // Convert angular extents to rotation angles: φ′ = φᶜ ± Δm/ζ
+            double phi_plus = phi_c + delta_m / zeta;
+            double phi_minus = phi_c - delta_m / zeta;
+
+            // Convert phi angles from radians to degrees before using scan parameters
+            double phi_plus_deg = phi_plus * 180.0 / M_PI;
+            double phi_minus_deg = phi_minus * 180.0 / M_PI;
+
+            // Transform rotation angles to image numbers using scan parameters
+            double z_plus =
+              image_range_start - 1 + ((phi_plus_deg - osc_start) / osc_width);
+            double z_minus =
+              image_range_start - 1 + ((phi_minus_deg - osc_start) / osc_width);
+
+            // Clamp to the actual image range and use floor/ceil for integer bounds
+            bbox.z_min =
+              std::max(image_range_start - 1,
+                       static_cast<int>(std::floor(std::min(z_plus, z_minus))));
+            bbox.z_max = std::min(
+              image_range_end, static_cast<int>(std::ceil(std::max(z_plus, z_minus))));
+        } else {
+            // Handle degenerate case where reflection is parallel to rotation axis
+            // In this case, the reflection spans the entire image range
+            bbox.z_min = image_range_start;
+            bbox.z_max = image_range_end;
+        }
+
+        extents.push_back(bbox);
+    }
+
+    return extents;
+}