Keep the constants out of the dual numbers, and clear only the boxes

Two costs in the per-image loop, each measured before it was touched.

Geometry refinement is the largest item in that loop - about half to two thirds of
its processor time on the datasets where the loop matters, and all of it on the host.
It runs three solves per image, and each one spends four fifths of itself inside the
solver at barely two iterations: the cost is not convergence, it is what every
residual evaluation does. The residual carried the blocks it does not refine as dual
numbers, so each evaluation recomputed the two detector rotations, the whole
orthogonalisation matrix, three cross products and the cell volume - all of them
constant for the image - through the derivative machinery, several million times per
run. Split the observed and predicted sides so the un-refined blocks pass as plain
doubles, evaluate the cell side once when the functor is built, and let the rotator
take a point whose type differs from the angle's. A dual number times a double is a
dual number times a dual number whose derivatives are zero, so the arithmetic is the
same one with the zeros removed.

Integration cleared the owner and mask images for the whole frame before every image.
On a large detector that is more than three hundred megabytes of writes to reset
pixels of which about one in twenty-five is ever marked, and it cost most of what the
integration kernels themselves cost. The marking kernel gained an unmarking mode - one
kernel, so the two cannot drift apart - and the engine clears whichever way is cheaper
for the frame in front of it, with a flag to force the full clear the first time and
after anything threw. The size test is not decoration: without it, clearing box by box
is slower than the memset on a small detector with many predictions, which is what the
measurement said before it was added.

Faster on thirteen of thirteen matched pairs: refinement by a quarter to a third,
whole-run wall by one to eight per cent depending on how much of the run is the loop.
The two changes pay in opposite regimes - refinement where the loop is processor-bound,
the clear where the detector is large enough for the card to be the constraint.

Every reflection file over seven datasets is byte-identical, and the solver did not
merely land in the same place: it took the same path, agreeing digit for digit on
iteration, residual and Jacobian evaluation counts. A new test runs two mismatched
frames through one engine and compares against a fresh one, which is what a mark left
behind would break.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016NNnL26LAvruQ9eLUUWvrJ

tests/BraggIntegrationEngineGPUTest.cpp

@@ -305,6 +305,61 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") {
     }
 }
 
+// The mask and the owner map are cleared by the run that marked them rather than at the start of the
+// next one, so a reused engine has to give the same answer as a fresh one. A first frame whose spots
+// are somewhere else entirely is what would show a leftover mark: a stale mask pixel is read as a
+// neighbour's signal and dropped from the background ring, a stale owner steals a pixel outright.
+TEST_CASE("BraggIntegrationEngineGPU_ReusedEngineMatchesFresh") {
+    if (get_gpu_count() == 0) {
+        WARN("No CUDA GPU present. Skipping BraggIntegrationEngineGPU_ReusedEngineMatchesFresh");
+        return;
+    }
+
+    for (OverlapMode ovl : {OverlapMode::Off, OverlapMode::Exclude}) {
+        const DiffractionExperiment experiment =
+                MakeExperiment(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, DetJF(2),
+                               0.0f, 0.0f, 0.0f, 0.0f, ovl);
+        const size_t width = experiment.GetXPixelsNum();
+        const size_t height = experiment.GetYPixelsNum();
+        const size_t npixel = experiment.GetPixelsNum();
+
+        // Two frames whose spot grids do not line up, so the first frame's marks fall on the second
+        // frame's background rings rather than back onto its own disks.
+        const Scene first = BuildScene(width, height, 47);
+        const Scene second = BuildScene(width, height, 60);
+        REQUIRE(first.predicted.size() > 60);
+        REQUIRE(second.predicted.size() > 60);
+
+        auto integrate = [&](BraggIntegrationEngineGPU &engine, const Scene &scene,
+                             const std::shared_ptr<CudaStream> &stream) {
+            ImagePreprocessorBufferGPU img(npixel);
+            for (size_t i = 0; i < npixel; ++i) img[i] = scene.image[i];
+            REQUIRE(cudaMemcpyAsync(img.getGPUBuffer(), img.getBuffer().data(),
+                                    npixel * sizeof(int32_t), cudaMemcpyHostToDevice, *stream) == cudaSuccess);
+            return engine.Run(img, scene.predicted, scene.predicted.size(), 7);
+        };
+
+        auto stream_fresh = std::make_shared<CudaStream>();
+        BraggIntegrationEngineGPU fresh(experiment, stream_fresh);
+        const auto out_fresh = integrate(fresh, second, stream_fresh);
+
+        auto stream_reused = std::make_shared<CudaStream>();
+        BraggIntegrationEngineGPU reused(experiment, stream_reused);
+        integrate(reused, first, stream_reused);
+        const auto out_reused = integrate(reused, second, stream_reused);
+
+        INFO("overlap mode " << static_cast<int>(ovl));
+        REQUIRE(out_reused.size() == out_fresh.size());
+        for (size_t i = 0; i < out_fresh.size(); ++i) {
+            INFO("reflection " << i);
+            CHECK(out_reused[i].h == out_fresh[i].h);
+            CHECK(out_reused[i].I == out_fresh[i].I);
+            CHECK(out_reused[i].sigma == out_fresh[i].sigma);
+            CHECK(out_reused[i].bkg == out_fresh[i].bkg);
+        }
+    }
+}
+
 // Hidden ([.]) benchmark: the raison d'etre of the GPU port is < 2 ms/frame (vs ~142 ms on the CPU
 // for ProfileIntegrate2D). Run explicitly with:  ./jfjoch_test "[bragg_bench]"
 TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") {