From 2c1123c267de2c31d40ad22b358f44232ab49ea9 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 28 Sep 2026 23:25:54 +0200 Subject: [PATCH] rugnux: calibration re-bin sums every image; clean out-of-host-memory failures Calibration (--mode calibration): RebinAndRefit gave each worker one AzimuthalIntegrationProfile that AzIntEngineCPU::Run clears on every image, so the re-binned profile held only the last image each worker read, added in completion order - on a 1800-frame LaB6 run the re-binned fit changed with -N (197 / 181 / 205 ring points at rms ~7 px). Each fixed block of images is now summed in image order into its own slot (ParallelBlocks) and the slots added in block order, so the profile is the sum over every image and the same at any thread count (rms 0.93 px, identical at -N 1, 7 and 32). On the LaB6 runs with a sane header the first pass still wins and the .poni is unchanged; with a header 100 px off the re-binned fit is now adopted (218 vs 140 ring points). New test Rugnux_CalibrationRebinsEveryImage spreads the rings over six frames and fails on the old code (241 vs 141 points at -N 1 vs -N 4). Host memory: out-of-memory now ends with "Processing failed: out of host memory (...) - this data set needs more RAM than is available" and exit 1. Tested with ulimit -v and an LD_PRELOAD allocator that refuses allocations from a chosen phase (CPU and GPU builds, image loop through merging and writing, calibration mode). Paths that crashed instead: - PostIndexingRefinement ran candidate blocks on bare std::threads, so a bad_alloc there called std::terminate (seen: "terminate called recursively", SIGABRT); now std::async futures. - FFTW aborts on its own failed allocation (CK(p) in kernel/alloc.c, also inside buffered transforms: BeamCenterFFTCPU, FFTIndexerCPU); rugnux now defines fftwf_assertion_failed to report it and exit 1. - libjpeg's default error_exit calls exit() from the diagnostic-JPEG thread ("Insufficient memory (case 12)"); WriteJPEGToMem now longjmps back and throws (bad_alloc for JERR_OUT_OF_MEMORY). - WorkerPool construction that fails to start a thread destroyed joinable threads (terminate); it now joins them and rethrows. Also: length_error counts as a fatal resource error, pinned host allocation failure is MemAllocFailed, and the GPU scaling fail-fast message says the CPU path needs a lot of host memory. Docs: very large cells and what running out of host memory looks like (including the OOM killer) in RUGNUX_INSTALL.md, pointer in RUGNUX.md. Validation: myob, cytc, lyso, sparse md5-identical to 82e6583e0 on both GPU and CPU builds. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- common/ParallelFor.h | 26 +++-- docs/CHANGELOG.md | 2 + docs/RUGNUX.md | 3 +- docs/RUGNUX_INSTALL.md | 35 ++++++ image_analysis/indexing/CUDAMemHelpers.h | 2 +- .../indexing/PostIndexingRefinement.cpp | 12 +- .../scale_merge/RotationScaleMergeGPU.cu | 6 +- preview/JFJochJPEG.cpp | 32 +++++- rugnux/Rugnux.cpp | 65 +++++------ rugnux/rugnux_cli.cpp | 44 +++++++- tests/JPEGTest.cpp | 7 ++ tests/RugnuxTest.cpp | 104 ++++++++++++++++++ 12 files changed, 280 insertions(+), 58 deletions(-) diff --git a/common/ParallelFor.h b/common/ParallelFor.h index 3f89857cd..276cbd35a 100644 --- a/common/ParallelFor.h +++ b/common/ParallelFor.h @@ -60,16 +60,26 @@ namespace parallel_detail { WorkerPool() { const unsigned hw = std::max(1u, std::thread::hardware_concurrency()); workers.reserve(hw - 1); - for (unsigned i = 0; i + 1 < hw; i++) // the submitting thread takes a share too - workers.emplace_back([this] { - // It inherited whatever CPUs the thread that first used the pool was kept to. - RestoreThreadAffinity(); - InWorkerFlag() = true; - Loop(); - }); + try { + for (unsigned i = 0; i + 1 < hw; i++) // the submitting thread takes a share too + workers.emplace_back([this] { + // It inherited whatever CPUs the thread that first used the pool was kept to. + RestoreThreadAffinity(); + InWorkerFlag() = true; + Loop(); + }); + } catch (...) { + // A thread that cannot be started (no memory left for its stack) fails the caller with + // that error - after joining the ones already running, since destroying a joinable + // std::thread calls std::terminate instead. + StopAndJoin(); + throw; + } } - ~WorkerPool() { + ~WorkerPool() { StopAndJoin(); } + + void StopAndJoin() { { std::lock_guard lock(m); stop = true; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 67b416ebf..a11d38730 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -16,6 +16,8 @@ * Rugnux's `--model` rigid-body refinement computes the model's structure factors from one copy of the model and derives the translation part of its Jacobian exactly, so it runs faster; placements can differ slightly. * `rugnux --model` also draws a null replicate again when its rigid-body placement ends near an orientation equivalent to the model's, not only when it starts near one. * Rugnux without a GPU no longer stops Bragg prediction at 20,000 reflections per frame, so a large unit cell integrates the same reflections on a CPU-only build as on a GPU. +* `rugnux --mode calibration` re-bins every image of a multi-image powder run for its second fit, not only the last image each thread read, and gives the same result at any thread count. +* Rugnux reports running out of host memory as `out of host memory ... - this data set needs more RAM than is available` and exits non-zero, instead of `std::bad_alloc` or a crash. ### 1.0.0-rc.173 diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index e8dc69f05..f2e6c5d82 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -110,7 +110,8 @@ A few things worth knowing before reaching for more flags: - **A run wants memory in proportion to what it integrates**, not to the detector: 2.5-14 GB of host RAM and 3-7 GB on the card over the datasets measured, both peaking in scaling and merging. [Installing Rugnux ▸ Memory](RUGNUX_INSTALL.md#memory) has the table and the two - flags that lower it. + flags that lower it. A very large cell needs far more, and a CPU-only build most of all — tens of + GB of host RAM ([Very large unit cells](RUGNUX_INSTALL.md#very-large-unit-cells)). - Everything else is in [Running Rugnux](RUGNUX_TUTORIAL.md#running-rugnux) and the full [Command-line options](RUGNUX_ADVANCED.md#command-line-options). diff --git a/docs/RUGNUX_INSTALL.md b/docs/RUGNUX_INSTALL.md index c99316c27..5daeaa01b 100644 --- a/docs/RUGNUX_INSTALL.md +++ b/docs/RUGNUX_INSTALL.md @@ -169,3 +169,38 @@ non-zero exit status. That is deliberate: a CUDA out-of-memory says nothing abou flight and everything about the machine, so skipping frames would leave a run that looks complete with quietly different merged numbers. A run is never silently short; it either completes or fails. On a GPU shared with another job, lower `-N` or wait for the card. + +### Very large unit cells + +The rule above holds for the cells in that table, and a very large cell sits far beyond it. A long +axis, a dense pattern and many frames multiply together: each frame predicts more reflections, each +reflection is split over more frames, and **tens of millions of partial observations are normal** for +a big cell. One sweep of a crystal with an axis of about 640 Å predicted 236 million reflections and +peaked at **17.5 GB of host memory on a GPU build and 72.5 GB on a CPU-only build**. A CPU-only build +keeps in host RAM what a GPU build keeps on the card, so it is the one that needs the large machine. +When such a set is too large for GPU scaling on the card, the run stops and says to run it on the +CPU — which then wants that much host memory. + +### When host memory runs out + +Nothing estimates the host memory a run will need before it starts, so a set that does not fit runs +until the memory is gone, and then ends in one of two ways. + +Where the allocation is refused — under an address-space limit (`ulimit -v`, which some batch +systems set), with strict overcommit, or for one request larger than the machine could ever grant — +the run stops with + +``` +Processing failed: out of host memory (std::bad_alloc) - this data set needs more RAM than is available +``` + +and a non-zero exit status. The CUDA driver reserves a large address space of its own, so a GPU build +under a tight `ulimit -v` can stop before reading anything, with `CUDA (GPU) error (out of memory)`. + +More often on Linux the allocation succeeds and the memory runs out as it is used. The kernel's +**OOM killer** then ends the process with SIGKILL, which no program can catch: the run stops with no +message of its own, the shell prints `Killed` (exit status 137), and `dmesg` or the batch system's +log has an `Out of memory: Killed process` line. A job capped by a cgroup memory limit ends the same +way. A run that disappears like this needs more RAM than the machine or the job was given, exactly as +if it had said so. `--no-merge` lowers the host peak, as above; otherwise the run needs a larger +machine. diff --git a/image_analysis/indexing/CUDAMemHelpers.h b/image_analysis/indexing/CUDAMemHelpers.h index 002a83f31..eeca489ea 100644 --- a/image_analysis/indexing/CUDAMemHelpers.h +++ b/image_analysis/indexing/CUDAMemHelpers.h @@ -303,7 +303,7 @@ public: CudaHostPtr() = default; explicit CudaHostPtr(size_t count) { if (cudaMallocHost(&ptr_, count * sizeof(T)) != cudaSuccess) - throw JFJochException(JFJochExceptionCategory::GPUCUDAError, + throw JFJochException(JFJochExceptionCategory::MemAllocFailed, "Failed to allocate pinned host memory"); } ~CudaHostPtr() { diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp index cd03b819d..d8d210eb0 100644 --- a/image_analysis/indexing/PostIndexingRefinement.cpp +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -6,7 +6,7 @@ #include #include -#include +#include namespace { struct config_ifssr final { @@ -160,15 +160,17 @@ std::vector Refine(const std::vector &in_spots, if (nblocks == 1) { RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr); } else { - std::vector workers; + // Futures, not bare threads: an exception in a block (out of memory) then reaches the caller + // instead of calling std::terminate. + std::vector> workers; workers.reserve(nblocks - 1); for (unsigned b = 1; b < nblocks; b++) - workers.emplace_back([&, b] { + workers.push_back(std::async(std::launch::async, [&, b] { RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr, b, nblocks); - }); + })); RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr, 0, nblocks); for (auto &w : workers) - w.join(); + w.get(); } std::vector candidates; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 38d52965a..14e19bdb7 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -654,12 +654,14 @@ namespace { size_t free_bytes = 0, total_bytes = 0; cudaMemGetInfo(&free_bytes, &total_bytes); cudaGetLastError(); - char msg[640]; + char msg[768]; snprintf(msg, sizeof(msg), "Scaling %lld partial observations needs more GPU memory than this card has: %.1f GB for " "their per-observation arrays alone (%zu bytes each), before the merge's other buffers, " "on a %.1f GB card. This data set is too large for GPU scaling on this card - run it on " - "the CPU (CUDA_VISIBLE_DEVICES= rugnux ..., or a CPU-only build)", + "the CPU (CUDA_VISIBLE_DEVICES= rugnux ..., or a CPU-only build) on a machine with plenty of " + "host memory: the CPU path keeps all of it in RAM, which for a set this large can be " + "tens of GB", static_cast(n_obs), double(n_obs) * OBS_BYTES / 1e9, OBS_BYTES, double(total_bytes) / 1e9); throw JFJochException(JFJochExceptionCategory::GPUCUDAError, msg); diff --git a/preview/JFJochJPEG.cpp b/preview/JFJochJPEG.cpp index 865cf98a1..26bc10fd4 100644 --- a/preview/JFJochJPEG.cpp +++ b/preview/JFJochJPEG.cpp @@ -2,15 +2,31 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#include #include #include #include +#include #include +#include #include "../common/JFJochException.h" #include "../common/CompressedImage.h" #include "JFJochJPEG.h" +namespace { + // libjpeg's own error handler prints its message and calls exit(), from whichever thread was + // encoding - out of memory included. This one jumps back to WriteJPEGToMem, which throws instead. + struct JPEGErrorManager { + jpeg_error_mgr mgr; + std::jmp_buf jump; + }; + + void JPEGErrorExit(j_common_ptr cinfo) { + std::longjmp(reinterpret_cast(cinfo->err)->jump, 1); + } +} + std::string WriteJPEGToMem(const CompressedImage& image, int quality) { unsigned char *buf = nullptr; @@ -24,9 +40,21 @@ std::string WriteJPEGToMem(const CompressedImage& image, int quality) { image.GetUncompressed(buffer); struct jpeg_compress_struct cinfo; - struct jpeg_error_mgr jerr; + JPEGErrorManager jerr; - cinfo.err = jpeg_std_error(&jerr); + cinfo.err = jpeg_std_error(&jerr.mgr); + jerr.mgr.error_exit = JPEGErrorExit; + if (setjmp(jerr.jump)) { + const bool out_of_memory = (jerr.mgr.msg_code == JERR_OUT_OF_MEMORY); + char message[JMSG_LENGTH_MAX]; + jerr.mgr.format_message(reinterpret_cast(&cinfo), message); + jpeg_destroy_compress(&cinfo); + free(buf); + if (out_of_memory) + throw std::bad_alloc(); + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + std::string("JPEG encoding failed: ") + message); + } jpeg_create_compress(&cinfo); jpeg_mem_dest(&cinfo, &buf, &buf_size); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 0fbb27295..a186f004b 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -147,7 +147,7 @@ namespace { // the dataset and still exits 0, so the run looks complete while every merged number downstream // has quietly changed. Fail the run instead, and let the operator free the card and repeat it. bool IsFatalResourceError(const std::exception &e) { - if (dynamic_cast(&e)) return true; + if (dynamic_cast(&e) || dynamic_cast(&e)) return true; const auto *jf = dynamic_cast(&e); return jf != nullptr && (jf->Category() == JFJochExceptionCategory::GPUCUDAError || jf->Category() == JFJochExceptionCategory::MemAllocFailed); @@ -747,42 +747,37 @@ std::optional Rugnux::RebinAndRefit(const CalibrationResult & .PoniRot2_rad(geometry.GetPoniRot2_rad()); const AzimuthalIntegrationMapping rebinned_mapping(rebinned_experiment, pixel_mask_); - std::atomic next_ordinal = 0; - std::mutex merge_m; - AzimuthalIntegrationProfile summed(rebinned_mapping); + // The sum over every image. The profile is a float sum, so the order it is added up in sets its + // last bits: each fixed block of images is folded in image order into its own slot, and the + // slots are added in block order, which gives the same profile at any thread count. + std::vector> block_sum(ReductionBlocks(images_to_process, 1)); + ParallelBlocks(images_to_process, std::max(1, config_.nthreads), [&](int b, int lo, int hi) { + std::vector decompression_buffer; + JFJochReaderRawImage img; + ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_); + ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum()); + AzIntEngineCPU azint(rebinned_mapping); + AzimuthalIntegrationProfile image_profile(rebinned_mapping); // cleared by every azint.Run + auto block = std::make_unique(rebinned_mapping); - const int nthreads = std::max(1, config_.nthreads); - std::vector> workers; - workers.reserve(nthreads); - for (int t = 0; t < nthreads; ++t) { - workers.emplace_back(std::async(std::launch::async, [&]() { - std::vector decompression_buffer; - JFJochReaderRawImage img; - ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_); - ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum()); - AzIntEngineCPU azint(rebinned_mapping); - AzimuthalIntegrationProfile local(rebinned_mapping); - - while (!cancelled_) { - const int ordinal = next_ordinal.fetch_add(1); - const int image_idx = start_image + ordinal * config_.stride; - if (image_idx >= end_image) break; - try { - if (!reader_.ReadRawImage(image_idx, img)) continue; - const uint8_t *image_ptr = img.image.GetUncompressedPtr(decompression_buffer); - preprocessor.Analyze(buffer, image_ptr, img.image.GetMode()); - azint.Run(buffer, local); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Error("Failed to re-integrate image {}: {}", image_idx, e.what()); - } + for (int ordinal = lo; ordinal < hi && !cancelled_; ++ordinal) { + const int image_idx = start_image + ordinal * config_.stride; + try { + if (!reader_.ReadRawImage(image_idx, img)) continue; + const uint8_t *image_ptr = img.image.GetUncompressedPtr(decompression_buffer); + preprocessor.Analyze(buffer, image_ptr, img.image.GetMode()); + azint.Run(buffer, image_profile); + *block += image_profile; + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Error("Failed to re-integrate image {}: {}", image_idx, e.what()); } - std::lock_guard lock(merge_m); - summed += local; - })); - } - for (auto &w : workers) - w.get(); + } + block_sum[b] = std::move(block); + }, 1); + AzimuthalIntegrationProfile summed(rebinned_mapping); + for (const auto &block : block_sum) + summed += *block; if (cancelled_) return std::nullopt; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 8eb64c648..71d662798 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -5,6 +5,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -13,7 +16,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -661,6 +666,37 @@ namespace { } return out; } + + // Running out of host memory ends a run with an exception that names only what failed in C++ terms - + // std::bad_alloc, a vector past its maximum size, a thread that could not get its stack - which + // does not tell the user what to do. The data set is not at fault; the machine is too small for it. + std::string FailureText(const std::exception &e) { + const auto *sys = dynamic_cast(&e); + const auto *jf = dynamic_cast(&e); + const bool out_of_memory = dynamic_cast(&e) != nullptr + || dynamic_cast(&e) != nullptr + || (sys != nullptr && sys->code() == std::errc::resource_unavailable_try_again) + || (jf != nullptr && jf->Category() == JFJochExceptionCategory::MemAllocFailed); + if (!out_of_memory) + return e.what(); + return fmt::format("out of host memory ({}) - this data set needs more RAM than is available", + e.what()); + } +} + +// FFTW ends the program itself when its allocator comes back empty (CK(p) in kernel/alloc.c) - while +// planning, or inside a buffered transform - from whichever thread that happens in, so it cannot become +// an exception. Defining its handler here keeps FFTW's own (kernel/assert.c) out of the link, and says +// what happened in the same words as every other way of running out of memory. +extern "C" void fftwf_assertion_failed(const char *s, int line, const char *file) { + std::fflush(stdout); + if (std::strcmp(s, "p") == 0) { + std::fprintf(stderr, "Processing failed: out of host memory (FFTW could not allocate) - " + "this data set needs more RAM than is available\n"); + std::_Exit(EXIT_FAILURE); + } + std::fprintf(stderr, "fftw: %s:%d: assertion failed: %s\n", file, line, s); + std::abort(); } // The body of main. Settings setters and the pipeline itself throw JFJochException on input the parser @@ -2151,7 +2187,7 @@ static int RunRugnux(int argc, char **argv) { try { result = process.Run(); } catch (const std::exception &e) { - logger.Error("Processing failed: {}", e.what()); + logger.Error("Processing failed: {}", FailureText(e)); exit(EXIT_FAILURE); } g_active_process = nullptr; @@ -2244,7 +2280,7 @@ static int RunRugnux(int argc, char **argv) { try { result = process.Run(); } catch (const std::exception &e) { - logger.Error("Calibration failed: {}", e.what()); + logger.Error("Calibration failed: {}", FailureText(e)); return 1; } g_active_process = nullptr; @@ -2812,7 +2848,7 @@ static int RunRugnux(int argc, char **argv) { WriteResultReport(output_prefix, input_file, experiment, failed, logger, provenance); exit(EXIT_FAILURE); } catch (const std::exception &e) { - logger.Error("Processing failed: {}", e.what()); + logger.Error("Processing failed: {}", FailureText(e)); exit(EXIT_FAILURE); } g_active_process = nullptr; @@ -2944,7 +2980,7 @@ int main(int argc, char **argv) { enable_gpu_numa_binding(); return RunRugnux(argc, argv); } catch (const std::exception &e) { - Logger("rugnux").Error("{}", e.what()); + Logger("rugnux").Error("{}", FailureText(e)); return EXIT_FAILURE; } } diff --git a/tests/JPEGTest.cpp b/tests/JPEGTest.cpp index f84a28e81..019eed6d9 100644 --- a/tests/JPEGTest.cpp +++ b/tests/JPEGTest.cpp @@ -33,3 +33,10 @@ TEST_CASE("JPEGTest","[JPEG]") { std::ofstream f("x.jpeg", std::ios::binary); f.write(s.data(), s.size()); } + +// libjpeg's own error handler calls exit(); an image it refuses must come back as an exception. +TEST_CASE("JPEGTest_ErrorThrows","[JPEG]") { + const std::vector val(70000); // wider than JPEG_MAX_DIMENSION + const CompressedImage image(val, 70000, 1); + CHECK_THROWS_AS(WriteJPEGToMem(image), JFJochException); +} diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 5eefe879b..70438d9e6 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -12,6 +12,7 @@ #include "../rugnux/Rugnux.h" #include "../rugnux/RugnuxCommandLine.h" #include "../rugnux/SpotWidth.h" +#include "../image_analysis/geom_refinement/Calibrants.h" namespace { // Write a small VDS dataset of `n` flat images and return nothing (prefix_master.h5 + @@ -135,6 +136,109 @@ TEST_CASE("Rugnux_Cancel", "[HDF5][Full]") { REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } +// A calibration re-integrates its images binned about the geometry it fitted, and fits again. That +// second profile must be the sum over every image: here each of six frames carries a different sixth +// of the LaB6 rings, so a profile missing a frame is missing rings. The header's beam centre is 5 px +// out, far enough that the first pass's smeared profile loses to the re-binned one - so the geometry +// that comes out is the one fitted on the re-binned sum, and it must not depend on the thread count. +TEST_CASE("Rugnux_CalibrationRebinsEveryImage", "[HDF5][DetGeomCalib]") { + RegisterHDF5Filter(); + constexpr int n = 6; + constexpr float true_beam_x = 1090.0f, true_beam_y = 1100.0f; + const std::vector rings_q = CalibrantRings("lab6"); + + DiffractionExperiment x(DetJF4M()); + x.FilePrefix("process_calib_in").ImagesPerTrigger(n).OverwriteExistingFiles(true); + x.BitDepthImage(16).ImagesPerFile(n).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true); + x.Compression(CompressionAlgorithm::NO_COMPRESSION); + x.BeamX_pxl(true_beam_x).BeamY_pxl(true_beam_y).DetectorDistance_mm(100).IncidentEnergy_keV(WVL_1A_IN_KEV) + .FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10)); + const DiffractionGeometry geom_true = x.GetDiffractionGeometry(); + x.BeamX_pxl(true_beam_x + 5.0f); // what the header says + + // Radii of the rings the beam reaches, in pixels about the true centre. + std::vector ring_radius; + for (const float q : rings_q) { + const float d = static_cast(2.0 * PI) / q; + if (d > geom_true.GetWavelength_A() / 2.0f) + ring_radius.push_back(geom_true.ResToPxl(d)); + } + + { + StartMessage start_message; + x.FillMessage(start_message); + FileWriter file_set(start_message); + ScanResultGenerator generator(x); + std::vector image(x.GetPixelsNum()); + for (int i = 0; i < n; i++) { + for (int64_t p = 0; p < x.GetPixelsNum(); p++) { + const float px = static_cast(p % x.GetXPixelsNum()); + const float py = static_cast(p / x.GetXPixelsNum()); + const float r = std::hypot(px - true_beam_x, py - true_beam_y); + float v = 5.0f; + for (size_t k = i; k < ring_radius.size(); k += n) { + const float dr = r - ring_radius[k]; + v += 400.0f * std::exp(-0.5f * dr * dr / (1.5f * 1.5f)); + } + image[p] = static_cast(std::lround(v)); + } + DataMessage message{}; + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.number = i; + REQUIRE_NOTHROW(file_set.WriteHDF5(message)); + generator.Add(message); + } + EndMessage end_message; + end_message.max_image_number = n; + generator.FillEndMessage(end_message); + file_set.WriteHDF5(end_message); + file_set.Finalize(); + } + + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("process_calib_in_master.h5")); + auto dataset = reader.GetDataset(); + REQUIRE(dataset); + DiffractionExperiment experiment(dataset->experiment); + auto azint = experiment.GetAzimuthalIntegrationSettings(); + azint.AzimuthalBinCount(32); + experiment.ImportAzimuthalIntegrationSettings(azint); + + std::vector results; + for (const int nthreads : {1, 4}) { + ProcessConfig config; + config.mode = ProcessMode::Calibration; + config.calibrant_ring_q = rings_q; + config.nthreads = nthreads; + config.write_process_h5 = false; + config.beam_center_check = false; + config.spot_finding = DiffractionExperiment::DefaultDataProcessingSettings(); + config.spot_finding.indexing = false; + + Rugnux process(reader, experiment, *dataset->pixel_mask, config); + ProcessResult result; + REQUIRE_NOTHROW(result = process.Run()); + REQUIRE(result.calibration.has_value()); + results.push_back(*result.calibration); + } + + for (const auto &cal : results) { + CHECK(cal.converged); + CHECK(cal.geometry.GetBeamX_pxl() == Catch::Approx(true_beam_x).margin(0.5)); + CHECK(cal.geometry.GetBeamY_pxl() == Catch::Approx(true_beam_y).margin(0.5)); + } + CHECK(results[0].ring_points == results[1].ring_points); + CHECK(results[0].rms_radial_pxl == results[1].rms_radial_pxl); + CHECK(results[0].geometry.GetBeamX_pxl() == results[1].geometry.GetBeamX_pxl()); + CHECK(results[0].geometry.GetBeamY_pxl() == results[1].geometry.GetBeamY_pxl()); + CHECK(results[0].geometry.GetDetectorDistance_mm() == results[1].geometry.GetDetectorDistance_mm()); + + reader.Close(); + remove("process_calib_in_master.h5"); + remove("process_calib_in_data_000001.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + TEST_CASE("RugnuxCommandLine_Full", "[process]") { DiffractionExperiment x(DetJF(1)); IndexingSettings idx;