diff --git a/common/DiffractionSpot.cpp b/common/DiffractionSpot.cpp index bb0c473f..ba58bbe2 100644 --- a/common/DiffractionSpot.cpp +++ b/common/DiffractionSpot.cpp @@ -7,8 +7,8 @@ DiffractionSpot::DiffractionSpot(uint32_t col, uint32_t line, int64_t in_photons) { if (in_photons < 0) in_photons = 0; - x = col * static_cast(in_photons); - y = line * static_cast(in_photons); + x = static_cast(col) * in_photons; + y = static_cast(line) * in_photons; pixel_count = 1; photons = in_photons; max_photons = in_photons; @@ -23,18 +23,10 @@ DiffractionSpot& DiffractionSpot::operator+=(const DiffractionSpot &other) { return *this; } -DiffractionSpot::DiffractionSpot(float x_sum, float y_sum, int64_t in_pixel_count, +DiffractionSpot::DiffractionSpot(int64_t x_sum, int64_t y_sum, int64_t in_pixel_count, int64_t in_photons, int64_t in_max_photons) : x(x_sum), y(y_sum), pixel_count(in_pixel_count), photons(in_photons), max_photons(in_max_photons) {} -DiffractionSpot::DiffractionSpot(const SpotToSave &save) { - x = save.x * static_cast(save.intensity); - y = save.y * static_cast(save.intensity); - pixel_count = 1; - photons = std::lround(save.intensity); - max_photons = save.maxc; -} - int64_t DiffractionSpot::Count() const { return photons; } @@ -46,7 +38,10 @@ int64_t DiffractionSpot::MaxCount() const { Coord DiffractionSpot::RawCoord() const { if (photons == 0) return {0, 0, 0}; - return {x / (float)photons, y / (float)photons, 0}; + // In double: the sums run past what a float mantissa holds on a bright spot, and the centroid is + // wanted to better than the last pixel bit. + return {static_cast(static_cast(x) / static_cast(photons)), + static_cast(static_cast(y) / static_cast(photons)), 0}; } int64_t DiffractionSpot::PixelCount() const { @@ -54,24 +49,38 @@ int64_t DiffractionSpot::PixelCount() const { } void DiffractionSpot::AddPixel(uint32_t col, uint32_t line, int64_t photons) { - this->x += col * (float) photons; - this->y += line * (float) photons; + this->x += static_cast(col) * photons; + this->y += static_cast(line) * photons; this->photons += photons; this->max_photons = std::max(this->max_photons, photons); this->pixel_count += 1; } void DiffractionSpot::ConvertToImageCoordinates(const DiffractionExperiment &experiment, uint16_t module_number) { - auto c_out = RawToConvertedCoordinate(experiment, module_number, RawCoord()); - this->x = c_out.x * (float) photons; - this->y = c_out.y * (float) photons; + // The raw -> image map is a signed axis swap (module axes are +/-1 along X or Y, never rotated - + // see DetectorGeometryModular::GetDirection) plus an integer translation: the module's origin in + // the assembled image, and the two-pixel multipixel gaps, which depend on where the centroid + // falls. Applied to the SUMS rather than to the centroid, it therefore leaves them exact + // integers; converting the centroid and multiplying it back by the photon count would put the + // rounding this class avoids straight back in. + const Coord centroid = RawCoord(); + const Coord fast = experiment.GetModuleFastDirection(module_number); + const Coord slow = experiment.GetModuleSlowDirection(module_number); + // Everything the map adds on top of the axis swap - origin plus gap corrections, an integer. + const Coord shift = RawToConvertedCoordinate(experiment, module_number, centroid) + - fast * centroid.x - slow * centroid.y; + + const int64_t x_raw = x, y_raw = y; + x = std::lround(shift.x) * photons + std::lround(fast.x) * x_raw + std::lround(slow.x) * y_raw; + y = std::lround(shift.y) * photons + std::lround(fast.y) * x_raw + std::lround(slow.y) * y_raw; } std::optional DiffractionSpot::Export(const DiffractionGeometry &geometry, int64_t image_num) const { if (photons == 0) return std::nullopt; - auto d = geometry.PxlToRes(x / (float) photons, y / (float) photons); + const Coord centroid = RawCoord(); + auto d = geometry.PxlToRes(centroid.x, centroid.y); float phi = 0.0f; if (geometry.GetRotation()) { @@ -81,8 +90,8 @@ std::optional DiffractionSpot::Export(const DiffractionGeometry &geo } return SpotToSave{ - .x = x / static_cast(photons), - .y = y / static_cast(photons), + .x = centroid.x, + .y = centroid.y, .phi = phi, .intensity = static_cast(photons), .maxc = max_photons, diff --git a/common/DiffractionSpot.h b/common/DiffractionSpot.h index e98fc34b..482f988c 100644 --- a/common/DiffractionSpot.h +++ b/common/DiffractionSpot.h @@ -10,19 +10,23 @@ // Definition of Bragg spot class DiffractionSpot { - float x = 0; - float y = 0; + // Photon-weighted position sums, sum(col * photons) and sum(line * photons). Integers, not a + // centroid: column, line and the per-pixel count are all integral, so the sums are exact and + // every implementation that builds them - host, GPU, whatever the build lets the compiler + // contract - arrives at the same bits. In float they did not: gcc fuses the multiply-add under + // -march=x86-64-v3 and not at the baseline, which left the last bit of the centroid a property + // of the build flags rather than of the data. + int64_t x = 0; + int64_t y = 0; int64_t pixel_count = 0; int64_t photons = 0; // total photon count int64_t max_photons = INT64_MIN; // maximum number of counts per pixel in the spot public: DiffractionSpot() = default; DiffractionSpot(uint32_t col, uint32_t line, int64_t photons); - DiffractionSpot(const SpotToSave &save); - // From already-summed quantities. x_sum/y_sum are sum(col * photons) / sum(line * photons), i.e. - // the members as AddPixel leaves them, NOT a centroid. Used by the GPU spot extractor, which - // computes those sums on the device. - DiffractionSpot(float x_sum, float y_sum, int64_t pixel_count, int64_t photons, int64_t max_photons); + // From already-summed quantities, i.e. the members as AddPixel leaves them, NOT a centroid. + // Used by the GPU spot extractor, which computes those sums on the device. + DiffractionSpot(int64_t x_sum, int64_t y_sum, int64_t pixel_count, int64_t photons, int64_t max_photons); DiffractionSpot& operator+=(const DiffractionSpot& spot); int64_t PixelCount() const; int64_t Count() const; diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 09a1ac7f..5c8d7c99 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -212,7 +212,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 for (int i = t; i < n; i += nthreads) atomicAdd(&count[label[root[i]]], 1); __syncthreads(); - // 3) sums, one thread per component, walking its members in ascending list order so the float + // 3) sums, one thread per component, walking its members in ascending list order so the // accumulation matches DiffractionSpot::AddPixel term for term. A component bigger than // max-pix is thrown away below, so it is not summed - which is also what keeps a whole lit // module or diffraction ring from turning into one thread walking tens of thousands of @@ -223,20 +223,16 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 const int want = count[l]; scratch[l].pixel_count = want; if (want > max_pix) continue; - float x = 0.0f, y = 0.0f; + long long x = 0, y = 0; long long photons = 0, max_photons = LLONG_MIN; int found = 0; for (int j = i; j < n && found < want; j++) { if (root[j] != static_cast(i)) continue; const long long counts = value[j]; - // Spelled out rather than left as "x += col * counts", because the rounding has to match - // DiffractionSpot::AddPixel term for term and the two compilers do not contract alike: - // gcc fuses AddPixel into a vfmadd whenever the build enables FMA (the CI flags do), so - // __fmaf_rn is the counterpart there. Where the host cannot fuse (a baseline -march, or - // MSVC, which does not contract by default) this leaves the last bit of the centroid - // differing by an ulp - see SpotExtractorGPUParityTest.cpp. - x = __fmaf_rn(static_cast(index[j] % width), static_cast(counts), x); - y = __fmaf_rn(static_cast(index[j] / width), static_cast(counts), y); + // Integers, exactly as DiffractionSpot::AddPixel does them, so host and device agree by + // construction - no rounding mode to match and nothing for either compiler to contract. + x += static_cast(index[j] % width) * counts; + y += static_cast(index[j] / width) * counts; photons += counts; max_photons = max(max_photons, counts); found++; diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index b3eafa2e..e00a86cc 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -21,9 +21,9 @@ // * both make a component's root its lowest list index, so both find the same roots; // * labels are handed out by a prefix sum over the roots in ascending order, which is the order the // host's second scan hands them out in, so the SPOT ORDER is identical; -// * the centroid sums are accumulated per component in ascending list order, in float, term for -// term as DiffractionSpot::AddPixel does them, with the rounding spelled out (see the comment at -// the sum itself - the two compilers do not contract a multiply-add alike). +// * the centroid sums are accumulated per component in ascending list order, in integers, term for +// term as DiffractionSpot::AddPixel does them, so there is no rounding for the two compilers to +// disagree about. // tests/SpotExtractorGPUParityTest.cpp holds the two to each other on realistic, occupancy-swept and // pathological frames, and checks that repeating a frame gives byte-identical output. @@ -38,8 +38,8 @@ // Per-component sums, in exactly the form DiffractionSpot holds them: x and y are sum(col*photons) // and sum(line*photons), not a centroid. struct SpotExtractorGPUSpot { - float x; - float y; + int64_t x; + int64_t y; int64_t photons; int64_t max_photons; int32_t pixel_count; diff --git a/tests/SpotExtractorGPUParityTest.cpp b/tests/SpotExtractorGPUParityTest.cpp index 622d30e4..b5a86d9b 100644 --- a/tests/SpotExtractorGPUParityTest.cpp +++ b/tests/SpotExtractorGPUParityTest.cpp @@ -195,26 +195,11 @@ public: } }; -// How many representable floats apart two values are. -int64_t UlpDistance(float a, float b) { - int32_t ia, ib; - memcpy(&ia, &a, sizeof(ia)); - memcpy(&ib, &b, sizeof(ib)); - if (ia < 0) ia = INT32_MIN - ia; // map to a monotone ordering across the sign - if (ib < 0) ib = INT32_MIN - ib; - return std::abs(static_cast(ia) - static_cast(ib)); -} - -// Everything that decides which spots exist and what they weigh is compared EXACTLY: the number of -// spots, their order, and each one's pixel count, photon sum and maximum. Those are integers, and a -// difference in any of them is a difference in the partition. -// -// The centroid is a float sum, and its last bit is a property of the BUILD rather than of either -// implementation: gcc contracts DiffractionSpot::AddPixel into an FMA whenever the flags allow it -// (the CI -march=x86-64-v3 does), while a baseline -march, or MSVC with its default /fp:precise, -// cannot. The extractor uses __fmaf_rn, so it is bit-exact against a host that fuses and one ulp -// off one that does not. An ulp bound catches a real divergence - which moves a centroid by pixels, -// not by 1e-4 of one - while staying true whichever way the host was built. +// Everything is compared EXACTLY: the number of spots, their order, and each one's pixel count, +// photon sum, maximum and centroid. All of them are built out of integer sums on both sides, so +// there is nothing here that a build flag can move - which is the point of accumulating in integers +// rather than in float, where gcc contracted the multiply-add under -march=x86-64-v3 and not at the +// baseline and left the last bit of the centroid a property of how the host was built. void RequireIdentical(const std::string &what, const std::vector &cpu, const std::vector &gpu) { @@ -228,13 +213,8 @@ void RequireIdentical(const std::string &what, // RawCoord divides the sums by the photon count, so comparing it compares the sums; a spot // whose photons sum to zero reports (0,0) on both sides by the same branch. const Coord c = cpu[i].RawCoord(), g = gpu[i].RawCoord(); -#ifdef __FMA__ REQUIRE(memcmp(&g.x, &c.x, sizeof(float)) == 0); REQUIRE(memcmp(&g.y, &c.y, sizeof(float)) == 0); -#else - REQUIRE(UlpDistance(g.x, c.x) <= 2); - REQUIRE(UlpDistance(g.y, c.y) <= 2); -#endif } }