mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

Rotation merge: drop rocking events with an overloaded pixel; capture uncertainty in the merge variance

Browse Source 
A saturated pixel in a spot means the brightest part of the reflection was
not measured. The integration used to drop the peak frame's partial (its
peak pixel is unreadable) and keep the flanks, so the combine extrapolated
the event from its tails by the partiality model: on a strongly
diffracting small-molecule crystal the strongest low-order reflections
read 2-3x low and were the largest SHELXL misfits. XDS drops such a
reflection (OVERLOAD); so does rugnux now.

- Integration (CPU + GPU engines): a reflection is `overloaded` when a
  signal-disk pixel is saturated, or unreadable on this frame but not in
  the run's pixel mask - EIGER/PILATUS write their error value for a
  pixel they could not count, which the preprocessor turns into a masked
  pixel like a gap's. The engines now receive the PixelMask to tell the
  two apart (an earlier attempt that re-classified the marker as
  saturation in the preprocessor broke a dataset whose gaps are not in
  the file's mask). An overloaded reflection is kept with its box sum,
  unfitted, only so its event can be recognised.
- Rotation combine (CPU + GPU): an event with any overloaded partial is
  dropped whole; counted in the log and the report
  (OBSERVATIONS_REJECTED_OVERLOAD=). The unmerged MTZ export drops it too.
- Everything else that reads reflections leaves an overloaded one out:
  AcceptReflection (stills merge, per-image scaling), the post-refinement
  gather, the axial-row sums.
- Capture uncertainty: the merge rebuilds each full's variance at the
  reflection's mean (counting_variance / ModelSigma) and dropped the
  capture term the combine had put into sigma, so a full extrapolated
  from part of its rocking curve merged at the weight of a whole one.
  Fulls now carry it (Obs::capture) and the rebuilt variance adds
  (capture * <I>)^2, host and device.

SHELXL R1 on rugnux's own integration (harness), median fix -> this:
citric acid .0648 -> .0420 (XDS .051; 221 events dropped, EXTI 1.02 -> 0.29),
HEPES .0396 -> .0381 (184), aspirin 20 keV .0387 -> .0385 (6),
aspirin 25 keV .0376 -> .0375 (5); metformin/nidppe/dnba/lalanine/cytidine
no overloads, unchanged. YAG .116 -> .128 (87 dropped; its scale loop does
not settle either way). Proteins and private subset: see the branch report.
Tests: BraggIntegrationEngineCPU_SaturatedPeakIsFlaggedNotDropped (new),
BraggIntegrationEngineGPU_MatchesCPU (overloaded flag compared),
AcceptReflection_ResolutionLimits, [write_reflections], [large].

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB
This commit is contained in:
leonarski_f
2026-10-04 21:01:40 +02:00
co-authored by Claude Opus 5.5
parent d061e52ba1
commit 37a8c8e24e
25 changed files with 252 additions and 80 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/BraggIntegrationEngineGPUTest.cpp
+15 -8
View File
@@ -2,6 +2,8 @@
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include <algorithm>
#include "../common/CUDAWrapper.h"
#ifdef JFJOCH_USE_CUDA
@@ -171,7 +173,7 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
ImagePreprocessorBuffer cpu_image(npixel);
for (size_t i = 0; i < npixel; ++i)
cpu_image[i] = scene.image[i];
BraggIntegrationEngineCPU cpu(experiment);
BraggIntegrationEngineCPU cpu(experiment, PixelMask(experiment));
const auto out_cpu = cpu.Run(cpu_image, scene.predicted, scene.predicted.size(), 5);
// GPU under test, identical input uploaded to the device
@@ -181,7 +183,7 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
gpu_image[i] = scene.image[i];
REQUIRE(cudaMemcpyAsync(gpu_image.getGPUBuffer(), gpu_image.getBuffer().data(),
npixel * sizeof(int32_t), cudaMemcpyHostToDevice, *stream) == cudaSuccess);
BraggIntegrationEngineGPU gpu(experiment, stream);
BraggIntegrationEngineGPU gpu(experiment, stream, PixelMask(experiment));
const auto out_gpu = gpu.Run(gpu_image, scene.predicted, scene.predicted.size(), 5);
// The ok/observed decisions are deterministic geometry, so both engines return the same set in
@@ -198,15 +200,20 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
ImagePreprocessorBuffer clean_image(npixel);
for (size_t i = 0; i < npixel; ++i)
clean_image[i] = clean_scene.image[i];
BraggIntegrationEngineCPU clean_cpu(experiment);
BraggIntegrationEngineCPU clean_cpu(experiment, PixelMask(experiment));
const auto out_clean = clean_cpu.Run(clean_image, clean_scene.predicted,
clean_scene.predicted.size(), 5);
CHECK(out_cpu.size() < out_clean.size());
// An unreadable pixel the (empty) mask does not explain reads as an overload, which keeps the
// reflection, flagged - so the cost shows in the measured ones.
const auto measured = std::count_if(out_cpu.begin(), out_cpu.end(),
[](const Reflection &r) { return !r.overloaded; });
CHECK(static_cast<size_t>(measured) < out_clean.size());
}
for (size_t i = 0; i < out_cpu.size(); ++i) {
INFO("mode " << static_cast<int>(mode) << " reflection " << i << " hkl " << out_cpu[i].h);
CHECK(out_gpu[i].h == out_cpu[i].h);
CHECK(out_gpu[i].image_number == out_cpu[i].image_number);
CHECK(out_gpu[i].overloaded == out_cpu[i].overloaded);
CHECK(out_gpu[i].bkg == Catch::Approx(out_cpu[i].bkg).epsilon(0.02).margin(0.5));
CHECK(out_gpu[i].I == Catch::Approx(out_cpu[i].I).epsilon(0.03).margin(2.0));
CHECK(out_gpu[i].sigma == Catch::Approx(out_cpu[i].sigma).epsilon(0.03).margin(0.5));
@@ -386,11 +393,11 @@ TEST_CASE("BraggIntegrationEngineGPU_ReusedEngineMatchesFresh") {
};
auto stream_fresh = std::make_shared<CudaStream>();
BraggIntegrationEngineGPU fresh(experiment, stream_fresh);
BraggIntegrationEngineGPU fresh(experiment, stream_fresh, PixelMask(experiment));
const auto out_fresh = integrate(fresh, second, stream_fresh);
auto stream_reused = std::make_shared<CudaStream>();
BraggIntegrationEngineGPU reused(experiment, stream_reused);
BraggIntegrationEngineGPU reused(experiment, stream_reused, PixelMask(experiment));
integrate(reused, first, stream_reused);
const auto out_reused = integrate(reused, second, stream_reused);
@@ -425,7 +432,7 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") {
REQUIRE(npixel == width * height);
auto stream = std::make_shared<CudaStream>();
BraggIntegrationEngineGPU gpu(experiment, stream);
BraggIntegrationEngineGPU gpu(experiment, stream, PixelMask(experiment));
for (int spacing : {28, 60}) {
const Scene scene = BuildScene(width, height, spacing);
const size_t nrefl = scene.predicted.size();
@@ -446,7 +453,7 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") {
const auto t1 = std::chrono::steady_clock::now();
const double ms = std::chrono::duration<double, std::milli>(t1 - t0).count() / iters;
BraggIntegrationEngineCPU cpu(experiment);
BraggIntegrationEngineCPU cpu(experiment, PixelMask(experiment));
ImagePreprocessorBuffer cpu_image(npixel);
for (size_t i = 0; i < npixel; ++i) cpu_image[i] = scene.image[i];
const auto c0 = std::chrono::steady_clock::now();