From 796b72f8b7f1e3c28d6f56dd85c0e4e4d6aa8a69 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 1 Nov 2026 08:56:31 +0100 Subject: [PATCH 1/4] rugnux: the header line and console summary print the wavelength --wavelength set The "Header:" log line and the console note under it printed the input file's wavelength (or the reader's 1 A placeholder where the file states none) even when --wavelength overrode it, so a run that processed, refined and wrote everything at the given wavelength announced another one at the top of its log. They now print the wavelength used, naming the file's value beside it when the option replaced it. --mode scale keeps the wavelength its reflections were integrated at, so it still prints the file's. The report (WAVELENGTH=, XDS.INP block), the merged and unmerged MTZ, the mmCIF and the _process.h5 were already written from the overridden experiment. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- rugnux/rugnux_cli.cpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index cc30c1bb5..9fb9bac46 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2163,15 +2163,22 @@ int RugnuxInputRun::PrepareRugnuxInput() { output_provenance.collection_date = dataset->arm_date; { // The header as the file states it, before anything is measured: what a question about a run - // is answered from first. + // is answered from first. The wavelength is the one used: --wavelength, where given, is not + // refined later, so the file's value would be printed for a run that never used it (--mode scale + // keeps the wavelength its reflections were integrated at). const auto &x = dataset->experiment; const auto gonio = x.GetGoniometer(); const auto [direct_x, direct_y] = x.GetDiffractionGeometry().GetDirectBeam_pxl(); - logger.Info("Header: {} images, detector {} {} x {} px of {:.4f} mm, wavelength {:.5f} A, distance " + const bool wavelength_given = wavelength_A && mode != RugnuxMode::Scale; + const float used_wavelength_A = wavelength_given ? *wavelength_A : x.GetWavelength_A(); + const std::string wavelength_source = + wavelength_given ? fmt::format(" (--wavelength; the file states {:.5f} A)", x.GetWavelength_A()) : ""; + logger.Info("Header: {} images, detector {} {} x {} px of {:.4f} mm, wavelength {:.5f} A{}, distance " "{:.3f} mm, beam {:.2f} {:.2f} px (direct beam {:.2f} {:.2f}), tilt {:.4f} {:.4f} deg, {}, " "polarization {}", reader.GetNumberOfImages(), x.GetDetectorDescription(), x.GetXPixelsNum(), x.GetYPixelsNum(), - x.GetPixelSize_mm(), x.GetWavelength_A(), x.GetDetectorDistance_mm(), x.GetBeamX_pxl(), + x.GetPixelSize_mm(), used_wavelength_A, wavelength_source, x.GetDetectorDistance_mm(), + x.GetBeamX_pxl(), x.GetBeamY_pxl(), direct_x, direct_y, x.GetPoniRot1_rad() * 180.0 / PI, x.GetPoniRot2_rad() * 180.0 / PI, gonio ? fmt::format("rotation from {:.3f} deg by {:.4f} deg per image", gonio->GetStart_deg(), @@ -2181,7 +2188,7 @@ int RugnuxInputRun::PrepareRugnuxInput() { : std::string("not stated")); console->OnNote(fmt::format("{} images, {} {} x {} px", reader.GetNumberOfImages(), x.GetDetectorDescription(), x.GetXPixelsNum(), x.GetYPixelsNum())); - console->OnNote(fmt::format("{:.4f} A, {:.1f} mm, {}", x.GetWavelength_A(), x.GetDetectorDistance_mm(), + console->OnNote(fmt::format("{:.4f} A, {:.1f} mm, {}", used_wavelength_A, x.GetDetectorDistance_mm(), gonio ? fmt::format("{:.3f} deg per image", gonio->GetIncrement_deg()) : std::string("no rotation axis"))); } From d67f10cda1c531df348bcede5777791a7262859d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 1 Nov 2026 08:56:31 +0100 Subject: [PATCH 2/4] rugnux: a _process.h5 carrying its images is typed as the images are stored Where the input is not a single HDF5 file the process file can link to - several masters read as one (JFJochMultiReader), CBF, CXI/XFEL input - the _process.h5 carries the images itself. Its start message was then typed by the experiment, which is pinned to signed 32-bit (the container images are DECODED into), while the images are passed on as stored: 16-bit signed from a JUNGFRAU in photons, unsigned 16/32-bit from an EIGER. The writer refuses an image that contradicts the declared format, so --write-process-h5 on two Cristallina (int16) runs failed with "Image is 16-bit signed, but the start message declares 32-bit signed". The pixel format, error and underload values are now taken from the first image as the reader hands it out, exactly as the linked-HDF5 case takes them from the source files. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- rugnux/Rugnux.cpp | 40 +++++++++++++++++++++------------------ tests/RugnuxTest.cpp | 45 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+), 18 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index cf9722375..d7fc2d31c 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1034,32 +1034,36 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // ask for the links only when the source is HDF5, and otherwise write the images into the file - // the per-image analysis is the point of it either way, and without pictures nothing downstream // could show the frame a reflection came from. + // Everything describing the pixels must describe the images the file links to or carries. experiment_ + // cannot: it is pinned to signed 32-bit, the container the reader DECODES images into, while the images + // themselves are passed on as stored - 16-bit signed from a JUNGFRAU in photons, unsigned 16- or 32-bit + // from an EIGER. Typed by experiment_, the virtual dataset would convert every value on read, and the + // images the file carries would contradict its header (refused by the writer). + StoredPixelFormat stored{static_cast(experiment_.GetByteDepthImage() * 8), experiment_.IsPixelSigned()}; if (auto *hdf5 = dynamic_cast(&reader_)) { // With the stride: image i of the process file is source image start_image + i * stride, so // the linked pictures line up with the per-image analysis written beside them. start_message.hdf5_source_data = hdf5->GetHDF5DataSource(start_image, images_to_process, config_.stride); - - // Everything describing the pixels must describe those linked files. experiment_ cannot: it - // is pinned to signed 32-bit, the container the reader hands images out in, which would type - // the virtual dataset int32 over unsigned 16- or 32-bit sources and silently convert every - // value on read. Only consumers going through the virtual view ever saw it - our own reader - // resolves the mapping and opens the source itself. - const auto stored = hdf5->GetStoredPixelFormat(); - start_message.bit_depth_image = stored.bit_depth; - // Unsigned only, for the reason FillMessage gives: DIALS remaps the top two codes whenever - // the field is present, and on a signed image those land inside the trusted range. - start_message.bit_depth_readout.reset(); - if (!stored.is_signed) - start_message.bit_depth_readout = stored.bit_depth; - start_message.pixel_signed = stored.is_signed; - start_message.error_value = stored.is_signed - ? -(int64_t(1) << (stored.bit_depth - 1)) // INTx_MIN - : (int64_t(1) << stored.bit_depth) - 1; // UINTx_MAX - start_message.underload_value = stored.is_signed ? start_message.error_value.value() + 1 : 0; + stored = hdf5->GetStoredPixelFormat(); } else { start_message.write_images = true; // no file to link to; the process file carries them + // As the reader hands the images out, which is how they are written. + JFJochReaderRawImage first; + if (reader_.ReadRawImage(start_image, first)) + stored = {static_cast(first.image.GetByteDepth() * 8), first.image.IsSigned()}; } + start_message.bit_depth_image = stored.bit_depth; + // Unsigned only, for the reason FillMessage gives: DIALS remaps the top two codes whenever + // the field is present, and on a signed image those land inside the trusted range. + start_message.bit_depth_readout.reset(); + if (!stored.is_signed) + start_message.bit_depth_readout = stored.bit_depth; + start_message.pixel_signed = stored.is_signed; + start_message.error_value = stored.is_signed + ? -(int64_t(1) << (stored.bit_depth - 1)) // INTx_MIN + : (int64_t(1) << stored.bit_depth) - 1; // UINTx_MAX + start_message.underload_value = stored.is_signed ? start_message.error_value.value() + 1 : 0; // The full-detector mask is 4 B/px and the ROI map 2 B/px - tens of megabytes each - and nothing // but the process file reads them. Copy them where the start message is handed on rather than here: diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 7114e22f3..dc257c7b9 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -12,6 +12,7 @@ #include "../common/ScanResultGenerator.h" #include "../writer/FileWriter.h" #include "../reader/JFJochHDF5Reader.h" +#include "../reader/JFJochMultiReader.h" #include "../rugnux/Rugnux.h" #include "../rugnux/RugnuxCommandLine.h" #include "../rugnux/RugnuxOutputFiles.h" @@ -94,6 +95,50 @@ TEST_CASE("Rugnux_AzInt", "[HDF5][Full]") { REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } +TEST_CASE("Rugnux_ProcessFileCarriesImagesAsStored", "[HDF5][Full]") { + // A reader with no HDF5 file to link to (here several inputs read as one) has its images written + // into the _process.h5, typed as they are stored - 16-bit signed - not as the 32-bit container + // they are decoded into. + WriteTestDataset("process_carry_in", 4); + + std::vector> inputs; + auto input = std::make_unique(); + REQUIRE_NOTHROW(input->ReadFile("process_carry_in_master.h5")); + inputs.push_back(std::move(input)); + JFJochMultiReader reader(std::move(inputs)); + auto dataset = reader.GetDataset(); + REQUIRE(dataset); + + ProcessConfig config; + config.mode = ProcessMode::AzimuthalIntegration; + config.nthreads = 2; + ProcessFileSink sink("process_carry_out"); + config.image_sink = &sink; + + Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config); + ProcessResult result; + REQUIRE_NOTHROW(result = process.Run()); + CHECK(result.images_processed == 4); + + { + JFJochHDF5Reader out; + REQUIRE_NOTHROW(out.ReadFile("process_carry_out_process.h5")); + CHECK(out.GetNumberOfImages() == 4); + CHECK(out.GetStoredPixelFormat().bit_depth == 16); + CHECK(out.GetStoredPixelFormat().is_signed); + std::shared_ptr img; + REQUIRE_NOTHROW(img = out.LoadImage(0)); + REQUIRE(img); + CHECK(img->Image()[100] == 5); + } + + reader.Close(); + remove("process_carry_in_master.h5"); + remove("process_carry_in_data_000001.h5"); + remove("process_carry_out_process.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + TEST_CASE("Rugnux_NoOutput", "[HDF5][Full]") { WriteTestDataset("process_noout_in", 6); From bb94e97f8e3ce45bfc78656e1d7dec1af8d4e871 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 1 Nov 2026 11:22:03 +0100 Subject: [PATCH 3/4] rugnux: a _process.h5 of an input it cannot link to holds no image data Where the input is anything but a single HDF5 master the process file can link to (several masters read as one, CBF, XTC, CXI, ...), the _process.h5 no longer copies the images into /entry/data/data: the copy made the file as large as the data. It now holds only the per-image analysis, and the earlier fix (d67f10cda), which typed that copy as the images are stored, goes with it - the int16 multi-master CrMX run that failed there writes its file without images. The reader takes such a file (no /entry/data/data, no linked data files) as one with /entry/detector/number images: the master's per-image analysis is read as for an integrated master and the reflections as always, so --mode scale reads it like any other, while LoadImage and ReadRawImage return no image instead of failing on the missing pixels (the viewer then shows none). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/HDF5.md | 4 +- reader/HDF5MetadataSource.cpp | 141 +++++++++++++++++++--------------- reader/JFJochHDF5Reader.cpp | 5 +- reader/JFJochHDF5Reader.h | 3 + rugnux/Rugnux.cpp | 52 ++++++------- tests/RugnuxTest.cpp | 34 ++++---- 6 files changed, 127 insertions(+), 112 deletions(-) diff --git a/docs/HDF5.md b/docs/HDF5.md index e6908eab6..aee90d853 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -118,7 +118,9 @@ of the raw images — without duplicating terabytes of data. This is a particularly FAIR-friendly artefact: it can be shared or archived alongside (or instead of) the raw data to convey what is in a dataset and how it was processed, while the `/entry/data/data` -VDS still resolves to the original images when they are available. `rugnux` can also process +VDS still resolves to the original images when they are available. Where the input is anything but a +single HDF5 master it could link to (several masters read as one, CBF, XTC, CXI, ...), the file has no +`/entry/data/data` at all - only the per-image analysis, which `--mode scale` reads all the same. `rugnux` can also process an equally-spaced *subset* of images (start/end/stride), producing a down-sampled reference set. ## 3. NXmx-standard content diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 0d6fa3b6d..f0c668375 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -538,6 +538,75 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen // Of a multichannel master; only its first channel is read std::string first_channel; + // The per-image analysis and metadata a master holds itself (VDS, integrated, or a _process.h5 that + // holds no images at all). + auto read_master_per_image = [&] { + if (master_file->Exists("/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image")) + dataset->efficiency = master_file->ReadVector( + "/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image"); + else + dataset->efficiency = std::vector(number_of_images, 1.0); + + if (master_file->Exists("/entry/roi")) + dataset->roi = master_file->FindLeafs("/entry/roi"); + + for (const auto &s: dataset->roi) { + dataset->roi_max.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/max")); + dataset->roi_sum.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/sum")); + dataset->roi_sum_sq.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/sum_sq")); + dataset->roi_npixel.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/npixel")); + dataset->roi_x.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/x")); + dataset->roi_y.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/y")); + } + + if (master_file->Exists("/entry/MX")) { + if (master_file->Exists("/entry/MX/peakCountUnfiltered")) + dataset->spot_count = master_file->ReadOptVector("/entry/MX/peakCountUnfiltered"); + else + dataset->spot_count = master_file->ReadOptVector("/entry/MX/nPeaks"); + + dataset->spot_count_low_res = master_file->ReadOptVector("/entry/MX/peakCountLowRes"); + dataset->spot_count_indexed = master_file->ReadOptVector("/entry/MX/peakCountIndexed"); + dataset->spot_count_ice_rings = master_file->ReadOptVector("/entry/MX/peakCountIceRingRes"); + dataset->spot_count_ice_control = master_file->ReadOptVector("/entry/MX/peakCountIceRingControl"); + + dataset->indexing_result = master_file->ReadOptVector("/entry/MX/imageIndexed"); + dataset->bkg_estimate = master_file->ReadOptVector("/entry/MX/bkgEstimate"); + dataset->spindle_blind_fraction = master_file->ReadOptVector("/entry/MX/spindleBlindFraction"); + dataset->ice_ring_score = master_file->ReadOptVector("/entry/MX/iceRingScore"); + dataset->resolution_estimate = master_file->ReadOptVector("/entry/MX/resolutionEstimate"); + dataset->profile_radius = master_file->ReadOptVector("/entry/MX/profileRadius"); + // Master files write indexedLatticeCount; data files / the per-file MX + // plugin use indexingLatticeCount. Accept either for backward compatibility. + dataset->indexing_lattice_count = master_file->ReadOptVector("/entry/MX/indexedLatticeCount"); + if (dataset->indexing_lattice_count.empty()) + dataset->indexing_lattice_count = master_file->ReadOptVector("/entry/MX/indexingLatticeCount"); + dataset->mosaicity_deg = master_file->ReadOptVector("/entry/MX/mosaicity"); + dataset->b_factor = master_file->ReadOptVector("/entry/MX/bFactor"); + dataset->image_scale_factor = master_file->ReadOptVector("/entry/MX/imageScaleFactor"); + dataset->image_scale_cc = master_file->ReadOptVector("/entry/MX/imageScaleCC"); + dataset->integrated_reflections = master_file->ReadOptVector("/entry/MX/integratedReflections"); + dataset->sweep_quality = master_file->ReadOptVector("/entry/MX/sweepQuality"); + if (master_file->Exists("/entry/MX/sweepQualityReasons")) { + const auto dim = master_file->GetDimension("/entry/MX/sweepQualityReasons"); + for (size_t i = 0; i < (dim.empty() ? 0 : dim[0]); i++) + dataset->sweep_quality_reasons.push_back( + master_file->ReadElement("/entry/MX/sweepQualityReasons", i) + .value_or("")); + } + dataset->frame_disposition = master_file->ReadOptVector("/entry/MX/frameDisposition"); + if (master_file->Exists("/entry/MX/frameDispositionCodes")) { + const auto dim = master_file->GetDimension("/entry/MX/frameDispositionCodes"); + for (size_t i = 0; i < (dim.empty() ? 0 : dim[0]); i++) + dataset->frame_disposition_codes.push_back( + master_file->ReadElement("/entry/MX/frameDispositionCodes", i) + .value_or("")); + } + } + if (master_file->Exists("/entry/image")) + dataset->max_value = master_file->ReadOptVector("/entry/image/max_value"); + }; + if (master_file->Exists("/entry/data/data")) { HDF5DataSet data_dataset(*master_file, "/entry/data/data"); HDF5Dcpl dcpl(data_dataset); @@ -554,70 +623,7 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen if (data_layout == HDF5DataSetLayout::VIRTUAL) vds_data_mappings = ReadVDSImageMappings(*master_file, "/entry/data/data"); - if (master_file->Exists("/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image")) - dataset->efficiency = master_file->ReadVector( - "/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image"); - else - dataset->efficiency = std::vector(number_of_images, 1.0); - - if (master_file->Exists("/entry/roi")) - dataset->roi = master_file->FindLeafs("/entry/roi"); - - for (const auto &s: dataset->roi) { - dataset->roi_max.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/max")); - dataset->roi_sum.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/sum")); - dataset->roi_sum_sq.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/sum_sq")); - dataset->roi_npixel.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/npixel")); - dataset->roi_x.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/x")); - dataset->roi_y.emplace_back(master_file->ReadVector("/entry/roi/" + s + "/y")); - } - - if (master_file->Exists("/entry/MX")) { - if (master_file->Exists("/entry/MX/peakCountUnfiltered")) - dataset->spot_count = master_file->ReadOptVector("/entry/MX/peakCountUnfiltered"); - else - dataset->spot_count = master_file->ReadOptVector("/entry/MX/nPeaks"); - - dataset->spot_count_low_res = master_file->ReadOptVector("/entry/MX/peakCountLowRes"); - dataset->spot_count_indexed = master_file->ReadOptVector("/entry/MX/peakCountIndexed"); - dataset->spot_count_ice_rings = master_file->ReadOptVector("/entry/MX/peakCountIceRingRes"); - dataset->spot_count_ice_control = master_file->ReadOptVector("/entry/MX/peakCountIceRingControl"); - - dataset->indexing_result = master_file->ReadOptVector("/entry/MX/imageIndexed"); - dataset->bkg_estimate = master_file->ReadOptVector("/entry/MX/bkgEstimate"); - dataset->spindle_blind_fraction = master_file->ReadOptVector("/entry/MX/spindleBlindFraction"); - dataset->ice_ring_score = master_file->ReadOptVector("/entry/MX/iceRingScore"); - dataset->resolution_estimate = master_file->ReadOptVector("/entry/MX/resolutionEstimate"); - dataset->profile_radius = master_file->ReadOptVector("/entry/MX/profileRadius"); - // Master files write indexedLatticeCount; data files / the per-file MX - // plugin use indexingLatticeCount. Accept either for backward compatibility. - dataset->indexing_lattice_count = master_file->ReadOptVector("/entry/MX/indexedLatticeCount"); - if (dataset->indexing_lattice_count.empty()) - dataset->indexing_lattice_count = master_file->ReadOptVector("/entry/MX/indexingLatticeCount"); - dataset->mosaicity_deg = master_file->ReadOptVector("/entry/MX/mosaicity"); - dataset->b_factor = master_file->ReadOptVector("/entry/MX/bFactor"); - dataset->image_scale_factor = master_file->ReadOptVector("/entry/MX/imageScaleFactor"); - dataset->image_scale_cc = master_file->ReadOptVector("/entry/MX/imageScaleCC"); - dataset->integrated_reflections = master_file->ReadOptVector("/entry/MX/integratedReflections"); - dataset->sweep_quality = master_file->ReadOptVector("/entry/MX/sweepQuality"); - if (master_file->Exists("/entry/MX/sweepQualityReasons")) { - const auto dim = master_file->GetDimension("/entry/MX/sweepQualityReasons"); - for (size_t i = 0; i < (dim.empty() ? 0 : dim[0]); i++) - dataset->sweep_quality_reasons.push_back( - master_file->ReadElement("/entry/MX/sweepQualityReasons", i) - .value_or("")); - } - dataset->frame_disposition = master_file->ReadOptVector("/entry/MX/frameDisposition"); - if (master_file->Exists("/entry/MX/frameDispositionCodes")) { - const auto dim = master_file->GetDimension("/entry/MX/frameDispositionCodes"); - for (size_t i = 0; i < (dim.empty() ? 0 : dim[0]); i++) - dataset->frame_disposition_codes.push_back( - master_file->ReadElement("/entry/MX/frameDispositionCodes", i) - .value_or("")); - } - } - if (master_file->Exists("/entry/image")) - dataset->max_value = master_file->ReadOptVector("/entry/image/max_value"); + read_master_per_image(); format = FileWriterFormat::NXmxVDS; } else if (master_file->Exists("/entry/data/data_000001")) { @@ -794,6 +800,13 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen image_size_x = master_file->GetInt("/entry/instrument/detector/detectorSpecific/x_pixels_in_detector"); image_size_y = master_file->GetInt("/entry/instrument/detector/detectorSpecific/y_pixels_in_detector"); number_of_images = 0; + // A _process.h5 of an input it could not link to holds the per-image analysis without the + // images: as many images as it has per-image entries, none of them with pixels. + if (master_file->Exists("/entry/detector/number")) { + number_of_images = master_file->GetDimension("/entry/detector/number")[0]; + images_per_file = number_of_images; + read_master_per_image(); + } } if (master_file->Exists("/entry/MX")) { diff --git a/reader/JFJochHDF5Reader.cpp b/reader/JFJochHDF5Reader.cpp index 4009a8cdf..22ffdb7ea 100644 --- a/reader/JFJochHDF5Reader.cpp +++ b/reader/JFJochHDF5Reader.cpp @@ -15,6 +15,7 @@ void JFJochHDF5Reader::ReadFile(const std::string &filename) { auto metadata = std::make_shared(); auto open_result = metadata->Open(filename, default_experiment); + has_pixels = open_result.image_layout.format != FileWriterFormat::NoFile; image_source_.Configure(std::move(open_result.image_layout)); // Original file: per-image metadata is co-located with the pixels. metadata->UseImageSourceForMetadata(&image_source_); @@ -67,6 +68,8 @@ bool JFJochHDF5Reader::ReadRawImage(int64_t image_number, JFJochReaderRawImage & if (!active_metadata_) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot load image if file not loaded"); + if (!has_pixels) + return false; auto loc = GetImageLocation(image_number); chunk = image_source_.PrepareDirectRead(loc); @@ -88,7 +91,7 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr &dataset std::unique_lock ul(hdf5_mutex); (void) update_dataset; - if (!dataset) + if (!dataset || !has_pixels) return false; if (!active_metadata_) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, diff --git a/reader/JFJochHDF5Reader.h b/reader/JFJochHDF5Reader.h index 175fedd89..053213873 100644 --- a/reader/JFJochHDF5Reader.h +++ b/reader/JFJochHDF5Reader.h @@ -31,6 +31,9 @@ class JFJochHDF5Reader : public JFJochReader { std::shared_ptr active_metadata_; std::string active_snapshot_; size_t number_of_images = 0; + // False for a file of per-image analysis without the images (a _process.h5 of an input it could + // not link to): its images have metadata, but no pixels to load. + bool has_pixels = false; bool LoadImage_i(std::shared_ptr &dataset, DataMessage &message, diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 482229d17..ae8213b00 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1035,41 +1035,35 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b start_message.file_format = FileWriterFormat::NXmxIntegrated; start_message.write_master_file = true; start_message.write_images = false; - // The process file LINKS to the original images rather than writing its own, which only an HDF5 - // input can offer: there is nothing for a virtual dataset to point at in a directory of CBFs. So - // ask for the links only when the source is HDF5, and otherwise write the images into the file - - // the per-image analysis is the point of it either way, and without pictures nothing downstream - // could show the frame a reflection came from. - // Everything describing the pixels must describe the images the file links to or carries. experiment_ - // cannot: it is pinned to signed 32-bit, the container the reader DECODES images into, while the images - // themselves are passed on as stored - 16-bit signed from a JUNGFRAU in photons, unsigned 16- or 32-bit - // from an EIGER. Typed by experiment_, the virtual dataset would convert every value on read, and the - // images the file carries would contradict its header (refused by the writer). - StoredPixelFormat stored{static_cast(experiment_.GetByteDepthImage() * 8), experiment_.IsPixelSigned()}; + // The process file LINKS to the original images rather than writing its own, which only a single + // HDF5 input can offer: there is nothing for a virtual dataset to point at in a directory of CBFs, + // an XTC stream or several masters read as one. Such a process file holds no /entry/data/data at + // all - a copy of the images would make it as large as the data, and the per-image analysis is + // the point of it. if (auto *hdf5 = dynamic_cast(&reader_)) { // With the stride: image i of the process file is source image start_image + i * stride, so // the linked pictures line up with the per-image analysis written beside them. start_message.hdf5_source_data = hdf5->GetHDF5DataSource(start_image, images_to_process, config_.stride); - stored = hdf5->GetStoredPixelFormat(); - } else { - start_message.write_images = true; // no file to link to; the process file carries them - // As the reader hands the images out, which is how they are written. - JFJochReaderRawImage first; - if (reader_.ReadRawImage(start_image, first)) - stored = {static_cast(first.image.GetByteDepth() * 8), first.image.IsSigned()}; + + // Everything describing the pixels must describe those linked files. experiment_ cannot: it + // is pinned to signed 32-bit, the container the reader hands images out in, which would type + // the virtual dataset int32 over unsigned 16- or 32-bit sources and silently convert every + // value on read. Only consumers going through the virtual view ever saw it - our own reader + // resolves the mapping and opens the source itself. + const auto stored = hdf5->GetStoredPixelFormat(); + start_message.bit_depth_image = stored.bit_depth; + // Unsigned only, for the reason FillMessage gives: DIALS remaps the top two codes whenever + // the field is present, and on a signed image those land inside the trusted range. + start_message.bit_depth_readout.reset(); + if (!stored.is_signed) + start_message.bit_depth_readout = stored.bit_depth; + start_message.pixel_signed = stored.is_signed; + start_message.error_value = stored.is_signed + ? -(int64_t(1) << (stored.bit_depth - 1)) // INTx_MIN + : (int64_t(1) << stored.bit_depth) - 1; // UINTx_MAX + start_message.underload_value = stored.is_signed ? start_message.error_value.value() + 1 : 0; } - start_message.bit_depth_image = stored.bit_depth; - // Unsigned only, for the reason FillMessage gives: DIALS remaps the top two codes whenever - // the field is present, and on a signed image those land inside the trusted range. - start_message.bit_depth_readout.reset(); - if (!stored.is_signed) - start_message.bit_depth_readout = stored.bit_depth; - start_message.pixel_signed = stored.is_signed; - start_message.error_value = stored.is_signed - ? -(int64_t(1) << (stored.bit_depth - 1)) // INTx_MIN - : (int64_t(1) << stored.bit_depth) - 1; // UINTx_MAX - start_message.underload_value = stored.is_signed ? start_message.error_value.value() + 1 : 0; // The full-detector mask is 4 B/px and the ROI map 2 B/px - tens of megabytes each - and nothing // but the process file reads them. Copy them where the start message is handed on rather than here: diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index dc257c7b9..fa971397c 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -95,15 +95,15 @@ TEST_CASE("Rugnux_AzInt", "[HDF5][Full]") { REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } -TEST_CASE("Rugnux_ProcessFileCarriesImagesAsStored", "[HDF5][Full]") { - // A reader with no HDF5 file to link to (here several inputs read as one) has its images written - // into the _process.h5, typed as they are stored - 16-bit signed - not as the 32-bit container - // they are decoded into. - WriteTestDataset("process_carry_in", 4); +TEST_CASE("Rugnux_ProcessFileWithoutImages", "[HDF5][Full]") { + // A reader with no single HDF5 file to link to (here several inputs read as one) gets a _process.h5 + // with the per-image analysis and no /entry/data/data, which reads back with its images' metadata + // and no pixels. + WriteTestDataset("process_noimg_in", 4); std::vector> inputs; auto input = std::make_unique(); - REQUIRE_NOTHROW(input->ReadFile("process_carry_in_master.h5")); + REQUIRE_NOTHROW(input->ReadFile("process_noimg_in_master.h5")); inputs.push_back(std::move(input)); JFJochMultiReader reader(std::move(inputs)); auto dataset = reader.GetDataset(); @@ -112,7 +112,7 @@ TEST_CASE("Rugnux_ProcessFileCarriesImagesAsStored", "[HDF5][Full]") { ProcessConfig config; config.mode = ProcessMode::AzimuthalIntegration; config.nthreads = 2; - ProcessFileSink sink("process_carry_out"); + ProcessFileSink sink("process_noimg_out"); config.image_sink = &sink; Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config); @@ -120,22 +120,22 @@ TEST_CASE("Rugnux_ProcessFileCarriesImagesAsStored", "[HDF5][Full]") { REQUIRE_NOTHROW(result = process.Run()); CHECK(result.images_processed == 4); + { + HDF5ReadOnlyFile file("process_noimg_out_process.h5"); + CHECK_FALSE(file.Exists("/entry/data/data")); + } { JFJochHDF5Reader out; - REQUIRE_NOTHROW(out.ReadFile("process_carry_out_process.h5")); + REQUIRE_NOTHROW(out.ReadFile("process_noimg_out_process.h5")); CHECK(out.GetNumberOfImages() == 4); - CHECK(out.GetStoredPixelFormat().bit_depth == 16); - CHECK(out.GetStoredPixelFormat().is_signed); - std::shared_ptr img; - REQUIRE_NOTHROW(img = out.LoadImage(0)); - REQUIRE(img); - CHECK(img->Image()[100] == 5); + CHECK(out.LoadImage(0) == nullptr); + CHECK_NOTHROW(out.ReadReflections(0, 3)); } reader.Close(); - remove("process_carry_in_master.h5"); - remove("process_carry_in_data_000001.h5"); - remove("process_carry_out_process.h5"); + remove("process_noimg_in_master.h5"); + remove("process_noimg_in_data_000001.h5"); + remove("process_noimg_out_process.h5"); REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } From 6e71f6fbe829e28b55cca73e4db3b085768447e8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 1 Nov 2026 11:31:53 +0100 Subject: [PATCH 4/4] rugnux: --mode scale on stills makes the merge of --mode mx, bit for bit The stills scale and merge (ScaleOnTheFly, partiality post-refinement, error model, merge, resolution cut, statistics) is one function, StillsScaleMerge, called by --mode mx and by --mode scale. On a stills _process.h5 the re-merge now reproduces the merged MTZ of the run that wrote it exactly (checked on a CrMX run with and without -S/-C, two CrMX masters read as one, an LCLS XTC set and a EuXFEL CXI set with an indexing ambiguity). Before, it differed by an overall scale of 0.87-0.89 and in the reflections kept. What it took: - Frames: the cell given with -C is stored (/entry/MX/stillsGivenUnitCell) and --mode scale leaves out the frames outside its population as --mode mx does (LeaveOutOtherCellForms, shared). With no reference MTZ it resolves the indexing ambiguity from the data, as --mode mx does before merging. - Reflections: excitation_error, observed, clipped and overloaded are written per reflection - the merge drops overloaded and weighs clipped reflections, and scale mode had neither - and lp/qe/flight are stored as double, so the reciprocal read back is the float the reflection held. A still keeps the resolution of its own lattice (no longer recomputed from the file's cell), its partiality and excitation error are no longer reset, and a mosaicity stored as NaN reads as none. - Rocking curve: the file now holds the curve the frames were predicted through, not the one the post-refinement balanced against the merge, which --mode scale then balanced a second time. - History: a space-group search merge leaves the frames' scales as it found them, so the merge that is written starts from the frames as integrated. RESULT-CHANGING for stills runs without -S (the final merge was seeded with the P1 search merge's scales): on a 3000-frame CrMX run CC1/2 0.9592 -> 0.9588, ISa 3.57 -> 3.58, the same reflections at an overall scale of 0.88. Rotation --mode scale is unchanged (to float rounding, from the exact correction factors). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- common/JFJochMessages.h | 4 +- docs/HDF5.md | 3 +- docs/RUGNUX_ADVANCED.md | 17 +- reader/HDF5MetadataSource.cpp | 38 ++- reader/JFJochReaderDataset.h | 4 +- rugnux/Rugnux.cpp | 2 +- rugnux/Rugnux.h | 32 ++- rugnux/RugnuxScaleMerge.cpp | 333 +++++++++++++----------- rugnux/rugnux_cli.cpp | 118 +++------ tests/JFJochReaderTest.cpp | 58 +++-- writer/HDF5DataFilePluginReflection.cpp | 20 +- writer/HDF5NXmx.cpp | 5 + 12 files changed, 363 insertions(+), 271 deletions(-) diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 390707bd1..0aa2ee4fb 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -473,8 +473,10 @@ struct EndMessage { // setting the merge works in - is only chosen afterwards, so the two settings are not always the // same one. Absent means they are, i.e. the identity. std::optional> reindex_matrix; - // Stills: the rocking curve the partialities were computed from (and the frames predicted through). + // Stills: the rocking curve the frames were predicted through (and the partialities computed from). std::optional stills_run_width; + // Stills: the cell given on the command line, whose population of frames the merge keeps. + std::optional stills_given_cell; // Geometry refined by the offline analysis. When present these overwrite the master-file geometry // at finalization (rugnux writes the refined values; the broker leaves them empty and keeps the diff --git a/docs/HDF5.md b/docs/HDF5.md index aee90d853..cca40f5f7 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -399,7 +399,8 @@ variants. | `reindexMatrix` | | change of basis from the setting the per-image data are in to the setting of `/entry/sample/unit_cell` (`[9]`, `int32`, flattened 3×3, row major) — see below | | `stillsRockingWidthFormula` | 1/Å, rad, relative | stills: the run's rocking width σ² = γ₀² + (γₑ d\*)² + (σ_bw λ d\*²/2)², as `[γ₀, γₑ, σ_bw]`; offline (Rugnux) only | | `stillsRadialScale` | | stills: the scale of the reciprocal lattice the rocking curve is centred on | -| `stillsRockingWidthTableDstar`, `stillsRockingWidthTable` | 1/Å | stills: the rocking width σ measured against d\*, which the partialities use (held beyond the first and last entry; the formula above is what the prediction window uses); `--mode scale` merges with formula and table instead of measuring another | +| `stillsRockingWidthTableDstar`, `stillsRockingWidthTable` | 1/Å | stills: the rocking width σ measured against d\* before the image pass, which the frames were predicted through (held beyond the first and last entry; the formula above is what the prediction window uses); the partiality post-refinement balances it against the merge, and `--mode scale` starts from this same curve instead of measuring another | +| `stillsGivenUnitCell` | Å, deg | stills: the cell given with `-C`, whose population of frames the merge kept (`[6]`); `--mode scale` leaves out the same frames | **Reindex matrix.** The per-image `h`, `k`, `l` and `latticeIndexed` are written as each image is processed, in the setting that image was *indexed* in. The space group is only settled afterwards, by diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index b32c70d43..39f53de61 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -108,10 +108,9 @@ How to resolve it: Two things to know about where the choice lands: -- **`--mode scale` cannot repair a stills run after the fact.** The per-image test happens at - integration time, so a `_process.h5` whose images were integrated without a reference MTZ has - already lost the distinction, and no re-merge brings it back; the data-only resolution - runs in the full pipeline, not in `--mode scale`. On rotation data, where +- **`--mode scale` resolves it as the full run does.** On stills with no reference MTZ it makes the + crystals consistent from the data alone before merging; the per-image test against a reference + happens at integration time and is kept in the stored reflections. On rotation data, where the ambiguity is one choice for the whole dataset, `--mode scale --model model.pdb` does resolve it. - **The choice reaches the written reflections**, not only the R-factors and the maps: the merged `.mtz` / `.cif` / `.hkl` (and `_unmerged.mtz`, where it is asked for) are written in the indexing @@ -291,9 +290,13 @@ declines it and the run stops with a message saying so. `--model` works here as that the model is placed by the rigid body only: a model it does not place, or a prediction, is not searched for by molecular replacement, `--mr` is refused, and the placed model is not refined. On rotation data, where the [indexing ambiguity](#the-indexing-ambiguity) is one choice for the whole dataset, it can settle the -enantiomorph label and the indexing — under the same fit test as always — and the re-merged reflections are written in what it settled. What no re-merge can -repair is a **stills** `_process.h5` integrated without a reference: there the ambiguity was -decided per image, at integration time, and the distinction is gone from the stored reflections. +enantiomorph label and the indexing — under the same fit test as always — and the re-merged reflections are written in what it settled. On +stills the re-merge is the full run's own merge: the same frames (a given cell's population, as +`/entry/MX/stillsGivenUnitCell` records it), the same rocking curve, the same data-only resolution of +the indexing ambiguity and the same scaling, so a +`_process.h5` merged again with the options it was written with reproduces the full run's reflections +(with a reference MTZ, the full run also lets the crystals' own cells overrule the reference's indexing +before merging, which `--mode scale` does not). Where the full run re-seated the lattice — the space group it settled on is in a different setting from the one each image was indexed in — the file records the change of basis as diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index f0c668375..096b95dfe 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -331,18 +331,24 @@ bool ReadReflectionsFromGroup(HDF5Object &file, auto bkg = file.ReadOptVector(image_group_name + "/background_mean"); // Written since the merge stopped back-deriving it; older _process.h5 do not carry it. auto var_bkg = file.ReadOptVector(image_group_name + "/background_variance"); - auto lp = file.ReadOptVector(image_group_name + "/lp"); + auto lp = file.ReadOptVector(image_group_name + "/lp"); // The sensor-efficiency part of lp, written since the correction existed. Absent in every // older file, and there the whole of lp is taken for Lorentz-polarization - which is what // those files mean, since nothing had corrected for the sensor when they were written. - auto qe = file.ReadOptVector(image_group_name + "/qe"); + auto qe = file.ReadOptVector(image_group_name + "/qe"); // The flight path, in the same reciprocal convention. Absent in every older file, and there // it is taken as 1: nothing had corrected for the flight path when those files were written. - auto flight = file.ReadOptVector(image_group_name + "/flight"); + auto flight = file.ReadOptVector(image_group_name + "/flight"); auto partiality = file.ReadOptVector(image_group_name + "/partiality"); auto phi = file.ReadOptVector(image_group_name + "/delta_phi"); auto zeta = file.ReadOptVector(image_group_name + "/zeta"); auto image_scale_corr = file.ReadOptVector(image_group_name + "/image_scale_corr"); + // Written since --mode scale started from the state the merge of --mode mx starts from; absent in + // older files, where they read as zero and false. + auto dist_ewald = file.ReadOptVector(image_group_name + "/excitation_error"); + auto observed = file.ReadOptVector(image_group_name + "/observed"); + auto clipped = file.ReadOptVector(image_group_name + "/clipped"); + auto overloaded = file.ReadOptVector(image_group_name + "/overloaded"); if (h.size() != l.size() || h.size() != k.size() || h.size() != d.size() || h.size() != predicted_x.size() || h.size() != predicted_y.size() @@ -361,16 +367,16 @@ bool ReadReflectionsFromGroup(HDF5Object &file, } float lp_val = 0.0; - if (lp.size() > i && lp[i] != 0.0f) - lp_val = 1.0f / lp[i]; + if (lp.size() > i && lp[i] != 0.0) + lp_val = static_cast(1.0 / lp[i]); float qe_val = 1.0f; - if (qe.size() > i && qe[i] != 0.0f) - qe_val = 1.0f / qe[i]; + if (qe.size() > i && qe[i] != 0.0) + qe_val = static_cast(1.0 / qe[i]); float flight_val = 1.0f; - if (flight.size() > i && flight[i] != 0.0f) - flight_val = 1.0f / flight[i]; + if (flight.size() > i && flight[i] != 0.0) + flight_val = static_cast(1.0 / flight[i]); float partiality_val = -1.0f; if (partiality.size() > i && partiality[i] >= 0.0f) @@ -418,12 +424,16 @@ bool ReadReflectionsFromGroup(HDF5Object &file, .bkg = bkg.at(i), .var_bkg = var_bkg_val, .sigma = int_err.at(i), + .dist_ewald = dist_ewald.size() > i ? dist_ewald[i] : 0.0f, .prescaling_corr = lp_val, .qe_corr = qe_val, .flight_corr = flight_val, .partiality = partiality_val, .zeta = zeta_val, - .image_scale_corr = image_scale_corr_val + .image_scale_corr = image_scale_corr_val, + .observed = observed.size() > i && observed[i], + .clipped = clipped.size() > i && clipped[i], + .overloaded = overloaded.size() > i && overloaded[i] }; reflections.emplace_back(r); } @@ -986,6 +996,8 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen } dataset->stills_run_width = w; } + if (const auto c = master_file->ReadOptVector("/entry/MX/stillsGivenUnitCell"); c.size() == 6) + dataset->stills_given_cell = UnitCell{.a = c[0], .b = c[1], .c = c[2], .alpha = c[3], .beta = c[4], .gamma = c[5]}; dataset->experiment.SampleName(master_file->GetString("/entry/sample/name")); @@ -1689,8 +1701,10 @@ std::vector HDF5MetadataSource::ReadReflections(size_t start // ── per-image mosaicity ─────────────────────────────────────────────── if (meta_file->Exists("/entry/MX/mosaicity")) { try { - outcome.mosaicity_deg = - meta_file->ReadElement("/entry/MX/mosaicity", meta_image_id); + // NaN is how the file says the image had none. + const auto mos = meta_file->ReadElement("/entry/MX/mosaicity", meta_image_id); + if (mos && std::isfinite(*mos)) + outcome.mosaicity_deg = mos; } catch (...) { } } diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index 53c898d33..b255e0ece 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -45,8 +45,10 @@ struct JFJochReaderDataset { // The e.s.d.s of the experiment's unit cell, a b c alpha beta gamma (/entry/sample/unit_cell_esd); empty // where the file has none. std::optional unit_cell_esd; - // Stills: the rocking curve the run's partialities were computed from (/entry/MX/stillsRockingWidth*). + // Stills: the rocking curve the run's frames were predicted through (/entry/MX/stillsRockingWidth*). std::optional stills_run_width; + // Stills: the cell the run was given, whose population of frames it merged (/entry/MX/stillsGivenUnitCell). + std::optional stills_given_cell; std::vector az_int_bin_to_q; std::vector az_int_bin_to_phi; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index ae8213b00..fbc18ffb2 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1311,7 +1311,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } if (config_.export_observations && !observation_verdicts_.empty()) result.outputs->observations = ObservationMtzTable(indexer->GetIntegrationOutcome(), observation_verdicts_, - observation_d_min_, end_msg.stills_run_width, + observation_d_min_, observation_run_width_, *result.consensus_cell, experiment_); if (config_.export_unmerged_partials) { if (observer) observer->OnPhase("Unmerged reflections"); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 8f0e88824..fc5eed035 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -356,6 +356,34 @@ bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const Vali std::pair ValidationEvidenceDifference(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate); +// Stills given a cell: only the frames whose refined cell belongs to its population (GivenCellPopulation) +// are merged - a second crystal form, fitted to the same cell, is left out, its reflections cleared. +// Returns the number of frames left out. +size_t LeaveOutOtherCellForms(std::vector &outcomes, const UnitCell &given_cell, Logger &logger); + +// The stills scale and merge, the one both --mode mx and --mode scale make: each image scaled against the +// merge of all, the partialities post-refined against it (unless for_search or --simple-stills), the error +// model fitted, the observations merged, cut at the resolution fall-off and the statistics taken. +// run_width is the rocking curve the frames were predicted through (measured here where there is none); +// verdicts asks for each observation's merge verdict (--export-observations). +struct StillsMergeResult { + std::vector merged; + MergeStatistics statistics; + std::optional reference_statistics; // --report-resolution + std::optional d_min; // the cut applied to merged + std::optional resolution_fit_A; + std::optional used_width; // the curve the post-refinement used + std::vector> verdicts; + double isa = 0.0, a = 0.0, b = 0.0; // error model, XDS convention + bool isa_resolved = true; +}; +StillsMergeResult StillsScaleMerge(const DiffractionExperiment &experiment, std::vector &outcomes, + const std::optional &reference_cell, + const std::optional &run_width, + const std::vector &reference_data, bool for_search, + bool verdicts, int nthreads, Logger &logger, + const std::function &phase); + // The longest real-space axis whose lattice rows this detector can still separate. Neighbouring rows of // an axis of length L lie 1/L apart in reciprocal space, which near the beam is D * lambda / L on the // detector; they stay apart while that exceeds the recorded spot width. So L = D * lambda / @@ -1613,9 +1641,11 @@ class Rugnux { // Stills: the rocking width measured before the image pass, on a sample of its frames, which every // frame is then predicted through (ImagePass) and the partialities are computed from (ScaleMerge). std::optional stills_run_width_; - // --export-observations: the written stills merge's verdict on every observation, and its resolution cut. + // --export-observations: the written stills merge's verdict on every observation, its resolution cut and + // the rocking curve its partialities were computed from. std::vector> observation_verdicts_; std::optional observation_d_min_; + std::optional observation_run_width_; static constexpr int STILLS_WINDOW_MIN_FRAMES = 2000; // the sample: a tenth of the run, at least this static constexpr int STILLS_WINDOW_FRACTION = 10; // Stills: where the frames' pixel values pile up, against the file's saturation limit (see there). diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 24f557e03..97ee3a3a9 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -599,6 +599,170 @@ static void WarnIfMetricHostsMoreSymmetry(const std::optional &cell, metric_system, metric_rotations, sg.xhm(), group_rotations); } +size_t LeaveOutOtherCellForms(std::vector &outcomes, const UnitCell &given_cell, Logger &logger) { + const double given = CrystalLattice(given_cell).CalcVolume(); + std::vector rel; + for (const auto &o : outcomes) + if (!o.reflections.empty()) + rel.push_back(o.latt.CalcVolume() / given - 1.0); + auto [low, high] = GivenCellPopulation(rel); + std::string all; + for (double r : rel) + all += fmt::format(" {:.4f}", r); + logger.Debug("Given cell: refined volumes relative to it:{}", all); + // A histogram valley alone does not make two crystal forms: one form whose cell drifts - a crystal + // swelling with dose, a cell refined a little differently in different orientations - falls into + // groups too (a drift through its range piles up at its two ends), and leaving one group out throws + // away part of the data of one structure. Measured: a rotation sweep processed as stills, its cell + // 0.25 % larger in volume in the second half than in the first, lost 841 of 1800 frames to the + // valley. Two forms of one protein that index on one cell differ by far more (9 % between the two + // forms of the battery's two-form sample), so groups closer than MIN_FORM_SEPARATION are one form. + // Their intensities would be the better judge, but a stills merge of groups taken from different + // orientations of one rotating crystal correlates only 0.8 with itself (0.66-0.88 measured), which + // is no clearer a line than the cells draw. + if (std::isfinite(low) || std::isfinite(high)) { + double sum_in = 0.0, sum_out = 0.0; + size_t n_in = 0, n_out = 0; + for (double r : rel) + (r >= low && r <= high) ? (sum_in += r, ++n_in) : (sum_out += r, ++n_out); + const double separation = std::fabs(sum_out / n_out - sum_in / n_in); + if (separation < MIN_FORM_SEPARATION) { + logger.Info("Given cell: the frames outside its volume population ({:+.2f}% to {:+.2f}%) are " + "{:.2f}% from it in volume, too close to be another crystal form - every frame is merged", + 100 * low, 100 * high, 100 * separation); + low = -std::numeric_limits::infinity(); + high = std::numeric_limits::infinity(); + } + } + size_t below = 0, above = 0; + for (auto &o : outcomes) { + if (o.reflections.empty()) + continue; + const double r = o.latt.CalcVolume() / given - 1.0; + if (r < low || r > high) { + o.reflections.clear(); + (r < low ? below : above)++; + } + } + if (below + above > 0) + logger.Info("Given cell: {} of {} indexed frames refined to a cell outside its population " + "({:+.1f}% to {:+.1f}% in volume) and are left out of the merge - {} smaller, {} larger", + below + above, rel.size(), 100 * low, 100 * high, below, above); + return below + above; +} + +StillsMergeResult StillsScaleMerge(const DiffractionExperiment &experiment, std::vector &outcomes, + const std::optional &reference_cell, + const std::optional &run_width, + const std::vector &reference_data, bool for_search, + bool verdicts, int nthreads, Logger &logger, + const std::function &phase) { + // Self-scale each image against the running merge with ScaleOnTheFly (fixed partiality), then merge + // directly. This runs even with an external reference: the reference is used only to break the + // per-image indexing ambiguity and to report CCref / inherit the R-free set, NOT as a scale anchor + // (scaling each image against a foreign dataset injects cross-dataset systematics and is a worse + // reference than the data's own merge). One pass: the per-image scale G is the exact one-pass + // solution, so iterating would only rebuild the reference from the freshly-scaled (noisy) data and + // degrade weak stills (measured CC1/2 collapse at the default 3 iters). + StillsMergeResult out; + // A search merge only looks: the frames leave it with the scales they came with, so the merge that is + // written starts from the frames as they were integrated - in --mode mx as in --mode scale, which reads + // them so from the _process.h5. + struct FrameScale { std::optional cc, g; std::optional cc_n; std::vector corr; }; + std::vector before; + if (for_search) + for (const auto &o : outcomes) { + FrameScale &f = before.emplace_back(FrameScale{o.image_scale_cc, o.image_scale_g, o.image_scale_cc_n, {}}); + for (const auto &r : o.reflections) + f.corr.push_back(r.image_scale_corr); + } + { + phase("Scaling images"); + TimingMark scale_mark(logger, "stills scaling images"); + ScaleOnTheFly(experiment, MergeAll(experiment, outcomes, nthreads, reference_cell)).Scale(outcomes, nthreads); + } + // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal + // orientation tilt against the running merge and recompute each reflection's partiality + scale + // correction (no re-integration). Never for the P1 search pass. + if (!for_search && experiment.GetScalingSettings().GetStillsPartialityRefine()) { + phase("Partiality post-refine"); + TimingMark refine_mark(logger, "stills partiality post-refine"); + StillsPartialityRefine refiner(experiment); + const double mean_tilt = refiner.Run(outcomes, nthreads, reference_cell, run_width, &out.used_width); + logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); + // The refine fits its OWN per-image scale, with no floor, and adopts it unconditionally + // when the image has no prior CC - which is exactly the state the guard above leaves. So + // the guard was protecting nothing: measured on a jet dataset, of 367 images it left + // unscaled only 8 were still unscaled in the output, and 53 reached the merge at or below + // a fiftieth of the run median, the worst at a 4525th. Whatever produced the final scale, + // it is the final scale that has to be credible. + ScaleOnTheFly::RejectCollapsedScales(outcomes); + } + + phase("Merging"); + TimingMark merge_mark(logger, "stills error model and merge"); + MergeOnTheFly merge_engine(experiment); + merge_engine.Threads(nthreads); + // For the de-novo P1 pass (for_search) drop the ice-ring reflections from the merged + // intensities and the error model, so the space-group search sees clean data; the final + // in-symmetry merge keeps them so completeness is not lost. + merge_engine.ExcludeIceRings(for_search); + if (reference_cell) + merge_engine.ReferenceCell(*reference_cell); + // Merge-consistency filter: on floods with many spurious lattices (e.g. XFEL large cells) + // most indexed crystals do not correlate with the true structure; --min-image-cc drops the + // crystals whose per-image CC to the reference is below the limit, so the merge keeps only + // the coherent (real) ones. Set before the error model so the model, the merge and the + // reported statistics are all fitted over the same set of images. The P1 search pass keeps + // every image: its job is to find the symmetry, not to produce final intensities. + merge_engine.FilterByImageCC(!for_search && experiment.GetScalingSettings().GetMinCCForImage() > 0.0); + merge_engine.RefineErrorModel(outcomes); + if (merge_engine.ErrorModelActive()) { + // XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is + // already the whole-range 1/b, there being no asymptote tier on this path. + const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} partiality c={:.3f} chi2={:.2f}", + em.a, em.b, em.isa, merge_engine.ErrorModelC(), merge_engine.ErrorModelChi2()); + } + merge_engine.AddImages(outcomes); + if (verdicts) + out.verdicts = merge_engine.Verdicts(outcomes); + out.merged = merge_engine.ExportReflections(); + merge_mark.Stop(); + + // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported + // shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full + // range, and the _process.h5 is written from the per-image outcomes, so no data is lost. A + // manual --scaling-high-resolution wins; the P1 search merge (for_search) is never cut so the + // space-group search still sees the full range. + const auto &cut_ss = experiment.GetScalingSettings(); + out.d_min = ApplyResolutionCutoff(out.merged, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(), + cut_ss.GetResolutionCCTarget(), for_search, logger, &out.resolution_fit_A); + + phase("Computing statistics"); + TimingMark stats_mark(logger, "stills statistics"); + out.statistics = merge_engine.MergeStats(out.merged, outcomes, reference_data, out.d_min); + // The same kept reflections binned over the reference range (--report-resolution); the + // stills error model is fitted over the whole range already, so no ISa of its own. + if (const auto range = cut_ss.GetReportResolutionRange_A(); range && !for_search) + out.reference_statistics = merge_engine.MergeStats(out.merged, outcomes, reference_data, std::nullopt, range); + // No asymptote tier, and the stills ISa already IS the whole-range 1/b. + const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); + out.isa = em.isa; + out.isa_resolved = merge_engine.ErrorModelBResolved(); + out.a = em.a; + out.b = em.b; + + for (size_t i = 0; i < before.size(); ++i) { + outcomes[i].image_scale_cc = before[i].cc; + outcomes[i].image_scale_g = before[i].g; + outcomes[i].image_scale_cc_n = before[i].cc_n; + for (size_t j = 0; j < before[i].corr.size(); ++j) + outcomes[i].reflections[j].image_scale_corr = before[i].corr[j]; + } + return out; +} + bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { auto &logger = p.logger; auto &result = p.result; @@ -695,57 +859,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { logger.Info("Sweep ends: left out {} frame(s) at the start and {} at the end that recorded no " "diffraction - a shutter or beam that was not on them", at_start, at_end); } - // Stills given a cell: only the frames whose refined cell belongs to its population are merged - // (see GivenCellPopulation) - a second crystal form, fitted to the same cell, is left out. + // Stills given a cell: only the frames of its population are merged. if (!experiment_.IsRotationIndexing() && stills_user_cell_) { - const double given = CrystalLattice(*stills_user_cell_).CalcVolume(); - std::vector rel; - for (const auto &o : outcomes) - if (!o.reflections.empty()) - rel.push_back(o.latt.CalcVolume() / given - 1.0); - auto [low, high] = GivenCellPopulation(rel); - std::string all; - for (double r : rel) - all += fmt::format(" {:.4f}", r); - logger.Debug("Given cell: refined volumes relative to it:{}", all); - // A histogram valley alone does not make two crystal forms: one form whose cell drifts - a crystal - // swelling with dose, a cell refined a little differently in different orientations - falls into - // groups too (a drift through its range piles up at its two ends), and leaving one group out throws - // away part of the data of one structure. Measured: a rotation sweep processed as stills, its cell - // 0.25 % larger in volume in the second half than in the first, lost 841 of 1800 frames to the - // valley. Two forms of one protein that index on one cell differ by far more (9 % between the two - // forms of the battery's two-form sample), so groups closer than MIN_FORM_SEPARATION are one form. - // Their intensities would be the better judge, but a stills merge of groups taken from different - // orientations of one rotating crystal correlates only 0.8 with itself (0.66-0.88 measured), which - // is no clearer a line than the cells draw. - if (std::isfinite(low) || std::isfinite(high)) { - double sum_in = 0.0, sum_out = 0.0; - size_t n_in = 0, n_out = 0; - for (double r : rel) - (r >= low && r <= high) ? (sum_in += r, ++n_in) : (sum_out += r, ++n_out); - const double separation = std::fabs(sum_out / n_out - sum_in / n_in); - if (separation < MIN_FORM_SEPARATION) { - logger.Info("Given cell: the frames outside its volume population ({:+.2f}% to {:+.2f}%) are " - "{:.2f}% from it in volume, too close to be another crystal form - every frame is merged", - 100 * low, 100 * high, 100 * separation); - low = -std::numeric_limits::infinity(); - high = std::numeric_limits::infinity(); - } - } - size_t below = 0, above = 0; - for (auto &o : outcomes) { - if (o.reflections.empty()) - continue; - const double r = o.latt.CalcVolume() / given - 1.0; - if (r < low || r > high) { - o.reflections.clear(); - (r < low ? below : above)++; - } - } - if (below + above > 0) { - logger.Info("Given cell: {} of {} indexed frames refined to a cell outside its population " - "({:+.1f}% to {:+.1f}% in volume) and are left out of the merge - {} smaller, {} larger", - below + above, rel.size(), 100 * low, 100 * high, below, above); + end_msg.stills_given_cell = stills_user_cell_; + if (LeaveOutOtherCellForms(outcomes, *stills_user_cell_, logger) > 0) { // The cell written is the mean over the frames merged, not over the other form's as well. result.consensus_cell = indexer->GetConsensusUnitCell(); end_msg.unit_cell = result.consensus_cell; @@ -1732,103 +1849,29 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { rsm->Run(for_search, /*full_stats=*/!geometry_prepass, /*measure_cc_before_corrections=*/measure_cc)); } - // Stills (rotation goes through RotationScaleMerge above): self-scale each image against the - // running merge with ScaleOnTheFly (fixed partiality), then merge directly. This runs even - // with an external reference: the reference is used only to break the per-image indexing - // ambiguity and to report CCref / inherit the R-free set, NOT as a scale anchor (scaling each - // image against a foreign dataset injects cross-dataset systematics and is a worse reference - // than the data's own merge). One pass: the per-image scale G is the exact one-pass solution, - // so iterating would only rebuild the reference from the freshly-scaled (noisy) data and - // degrade weak stills (measured CC1/2 collapse at the default 3 iters). - { - phase("Scaling images (" + label + ")"); - TimingMark scale_mark(logger, "stills scaling images " + label); - auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), config_.nthreads, - result.consensus_cell); - indexer->ScaleAllImages(merge_result); + // Stills (rotation goes through RotationScaleMerge above). + const bool verdicts = config_.export_observations && !for_search && write_output; + auto sm = StillsScaleMerge(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell, + stills_run_width_, config_.reference_data, for_search, verdicts, config_.nthreads, + logger, [&](const std::string &name) { phase(name + " (" + label + ")"); }); + if (verdicts) { + observation_verdicts_ = std::move(sm.verdicts); + observation_d_min_ = sm.d_min; + observation_run_width_ = sm.used_width; } - // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal - // orientation tilt against the running merge and recompute each reflection's partiality + scale - // correction (no re-integration). Never for the P1 search pass. - if (!for_search && experiment_.GetScalingSettings().GetStillsPartialityRefine()) { - phase("Partiality post-refine (" + label + ")"); - TimingMark refine_mark(logger, "stills partiality post-refine " + label); - StillsPartialityRefine refiner(experiment_); - const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads, result.consensus_cell, - stills_run_width_, &end_msg.stills_run_width); - logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); - // The refine fits its OWN per-image scale, with no floor, and adopts it unconditionally - // when the image has no prior CC - which is exactly the state the guard above leaves. So - // the guard was protecting nothing: measured on a jet dataset, of 367 images it left - // unscaled only 8 were still unscaled in the output, and 53 reached the merge at or below - // a fiftieth of the run median, the worst at a 4525th. Whatever produced the final scale, - // it is the final scale that has to be credible. - ScaleOnTheFly::RejectCollapsedScales(indexer->GetIntegrationOutcome()); - } - const std::vector &merge_input = indexer->GetIntegrationOutcome(); - - phase("Merging"); - TimingMark merge_mark(logger, "stills error model and merge " + label); - MergeOnTheFly merge_engine(experiment_); - merge_engine.Threads(config_.nthreads); - // For the de-novo P1 pass (for_search) drop the ice-ring reflections from the merged - // intensities and the error model, so the space-group search sees clean data; the final - // in-symmetry merge keeps them so completeness is not lost. - merge_engine.ExcludeIceRings(for_search); - if (result.consensus_cell.has_value()) - merge_engine.ReferenceCell(*result.consensus_cell); - // Merge-consistency filter: on floods with many spurious lattices (e.g. XFEL large cells) - // most indexed crystals do not correlate with the true structure; --min-image-cc drops the - // crystals whose per-image CC to the reference is below the limit, so the merge keeps only - // the coherent (real) ones. Set before the error model so the model, the merge and the - // reported statistics are all fitted over the same set of images. The P1 search pass keeps - // every image: its job is to find the symmetry, not to produce final intensities. - merge_engine.FilterByImageCC(!for_search - && experiment_.GetScalingSettings().GetMinCCForImage() > 0.0); - merge_engine.RefineErrorModel(merge_input); - if (merge_engine.ErrorModelActive()) - { - // XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is - // already the whole-range 1/b, there being no asymptote tier on this path. - const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); - logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} partiality c={:.3f} chi2={:.2f}", - em.a, em.b, em.isa, merge_engine.ErrorModelC(), merge_engine.ErrorModelChi2()); - } - merge_engine.AddImages(merge_input); - if (config_.export_observations && !for_search && write_output) - observation_verdicts_ = merge_engine.Verdicts(merge_input); - + // The curve the frames were predicted through, which --mode scale passes in again: the one the + // post-refinement used is balanced against this run's merge, and balancing it again would move it. + end_msg.stills_run_width = stills_run_width_; + if (sm.resolution_fit_A) + result.resolution_fit_A = sm.resolution_fit_A; + result.error_model_isa = sm.isa; + result.error_model_isa_resolved = sm.isa_resolved; + result.error_model_a = sm.a; + result.error_model_b = sm.b; ScaleMergeResult out; - out.merged = merge_engine.ExportReflections(); - merge_mark.Stop(); - - // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported - // shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full - // range, and the _process.h5 is written from the per-image outcomes, so no data is lost. A - // manual --scaling-high-resolution wins; the P1 search merge (for_search) is never cut so the - // space-group search still sees the full range. (Rotation is cut inside RotationScaleMerge.) - const auto &cut_ss = experiment_.GetScalingSettings(); - const std::optional effective_d_min = ApplyResolutionCutoff( - out.merged, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(), - cut_ss.GetResolutionCCTarget(), for_search, logger, &result.resolution_fit_A); - if (config_.export_observations && !for_search && write_output) - observation_d_min_ = effective_d_min; - - phase("Computing statistics"); - TimingMark stats_mark(logger, "stills statistics " + label); - out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min); - // The same kept reflections binned over the reference range (--report-resolution); the - // stills error model is fitted over the whole range already, so no ISa of its own. - if (const auto range = cut_ss.GetReportResolutionRange_A(); range && !for_search) - out.reference_statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, - std::nullopt, range); - { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b - const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); - result.error_model_isa = em.isa; - result.error_model_isa_resolved = merge_engine.ErrorModelBResolved(); - result.error_model_a = em.a; - result.error_model_b = em.b; - } + out.merged = std::move(sm.merged); + out.statistics = std::move(sm.statistics); + out.reference_statistics = std::move(sm.reference_statistics); out.cc_half_before_corrections = out.statistics.overall.cc_half; // no surfaces on this path logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label); Note(fmt::format("merge in {}: {} unique", label, out.merged.size())); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 0380dad44..421fb317b 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -65,6 +65,7 @@ #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/RfreeFlags.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" // MatchReference +#include "../image_analysis/scale_merge/ReferenceFreeReindex.h" #include "../image_analysis/scale_merge/ScaleOnTheFly.h" #include "../image_analysis/scale_merge/StillsPartialityRefine.h" #include "../image_analysis/scale_merge/SublatticeRelabel.h" @@ -2440,18 +2441,29 @@ int RugnuxInputRun::PrepareRugnuxInput() { } } + // Stills given a cell: the frames of its population, as --mode mx merged them. + if (!experiment.IsRotationIndexing() && dataset->stills_given_cell) + LeaveOutOtherCellForms(reflections, *dataset->stills_given_cell, logger); // Stills: a crystal indexed in a lattice containing the file's cell is put on that cell, as the full - // pipeline does - before the resolution of every reflection is recomputed from the cell below, - // which is right only for reflections indexed on it. + // pipeline does. if (!experiment.IsRotationIndexing()) { const auto relabel = RelabelOnSublattice(reflections, experiment.GetUnitCell().value(), experiment.GetIndexingSettings()); logger.Info("Put {} crystal(s) indexed in a lattice containing the file's cell on that cell", relabel.frames); } - auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections); - logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images); - experiment.ImagesPerTrigger(refl_stats.n_images); + // A rotation sweep's reflections are put on the one cell they are merged in; a still keeps the + // resolution of its own lattice, as --mode mx merges it. + if (experiment.IsRotationIndexing()) + UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections); + size_t n_reflections = 0, n_images = 0; + for (const auto &o : reflections) + if (!o.reflections.empty()) { + n_reflections += o.reflections.size(); + ++n_images; + } + logger.Info("Read {} reflections from {} images", n_reflections, n_images); + experiment.ImagesPerTrigger(n_images); // Ice-ring handling, as the full pipeline does it (Rugnux.cpp): flag reflections on a // hexagonal-ice powder ring so scaling skips them while the merge keeps them. The flag is not @@ -2516,7 +2528,7 @@ int RugnuxInputRun::PrepareRugnuxInput() { // Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine -> // merge, including the default-on decay + absorption correction surfaces. It does not support // external-reference scaling or wedge refinement. - // Everything else (stills, reference scaling) uses ScaleOnTheFly + MergeOnTheFly. + // Stills go through StillsScaleMerge, the merge of --mode mx. const bool is_rotation = experiment.IsRotationIndexing(); if (is_rotation) { if (!reference_data.empty() @@ -2549,80 +2561,30 @@ int RugnuxInputRun::PrepareRugnuxInput() { reference_statistics = std::move(r.reference_statistics); reference_isa = r.reference_isa; } else { - // Scaling self-references: the reference MTZ (if any) fixes the cell/space group, reports - // CCref and provides the R-free test set, but is NOT a scale anchor - scaling each image - // against a foreign dataset injects cross-dataset systematics and is a worse reference than - // the data's own merge. The per-image scale G is the exact one-pass solution, so one pass - // (iterating a self-rebuilt reference only re-fits the freshly-scaled noise on weak stills). - // The file holds each still as it was merged - with its post-refined partiality - and not the - // excitation errors its prediction made. Put every reflection back where the prediction left it, - // at p = 1 and with the excitation error of the stored lattice and geometry, so that the scaling - // and the post-refinement start from the state they start from in the full pipeline. - for (auto &o: reflections) { - const Coord S0 = o.geom.GetScatteringVector(); - const double inv_lambda = 1.0 / o.geom.GetWavelength_A(); - for (auto &r: o.reflections) { - const Coord S = o.latt.Astar() * static_cast(r.h) + o.latt.Bstar() * static_cast(r.k) - + o.latt.Cstar() * static_cast(r.l) + S0; - r.dist_ewald = static_cast(std::sqrt(S.x * S.x + S.y * S.y + S.z * S.z) - inv_lambda); - r.partiality = 1.0f; + // A group whose lattice admits alternative indexings, with no reference to settle them: the + // crystals were integrated each in one of them at random, and are made consistent from the + // data alone, as --mode mx makes them before its merge. + if (const auto sg = experiment.GetGemmiSpaceGroup(); sg && reference_data.empty()) { + const auto ops = ReindexAmbiguityOperators(*experiment.GetUnitCell(), *sg); + if (!ops.empty()) { + const auto rf = ResolveReindexReferenceFree(reflections, experiment, *experiment.GetUnitCell(), + ops, nthreads, /*apply=*/true); + logger.Info("Indexing ambiguity resolved from the data: {} of {} crystals reindexed", + rf.reindexed, rf.crystals); } } - ScaleOnTheFly(experiment, MergeAll(experiment, reflections, nthreads, experiment.GetUnitCell())).Scale(reflections, nthreads); - // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal - // orientation tilt against the merge and recompute each reflection's partiality + scale correction - // (no re-integration), then merge with the improved corrections. - if (experiment.GetScalingSettings().GetStillsPartialityRefine()) { - StillsPartialityRefine refiner(experiment); - // The curve the mx run computed its partialities from (and predicted the frames through), - // where the file carries it: measured again on frames predicted through it, it is another. - const double mean_tilt = refiner.Run(reflections, nthreads, experiment.GetUnitCell(), - dataset->stills_run_width); - logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); - // As in the full pipeline: the refine fits its own per-image scale, with no floor. - ScaleOnTheFly::RejectCollapsedScales(reflections); - } - MergeOnTheFly merge_engine(experiment); - merge_engine.Threads(nthreads); - merge_engine.ReferenceCell(experiment.GetUnitCell()); - // --min-image-cc has to hold for the merge itself, not only for the reported statistics. - merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0); - // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full - // pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity - // model and produce much worse stills intensities (no (b*I)^2 systematic term at all). - merge_engine.RefineErrorModel(reflections); - if (merge_engine.ErrorModelActive()) - { - // XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is - // already the whole-range 1/b, there being no asymptote tier on this path. - const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); - logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} partiality c={:.3f} chi2={:.2f}", - em.a, em.b, em.isa, merge_engine.ErrorModelC(), merge_engine.ErrorModelChi2()); - } - merge_engine.AddImages(reflections); - merged_reflections = merge_engine.ExportReflections(); - - // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual - // --scaling-high-resolution wins, otherwise trim the written reflections + reported shells - // to the CC1/2 fall-off. (Rotation is cut inside RotationScaleMerge above.) - const auto &cut_ss = experiment.GetScalingSettings(); - // The offline --scale path re-scales a stored _process.h5 and is never a P1 search merge. - const std::optional effective_d_min = ApplyResolutionCutoff( - merged_reflections, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(), - cut_ss.GetResolutionCCTarget(), /*for_search=*/false, logger, &resolution_fit_A); - - merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, - effective_d_min); - if (const auto range = cut_ss.GetReportResolutionRange_A()) - reference_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, - std::nullopt, range); - { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b - const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); - error_model_isa = em.isa; - error_model_isa_resolved = merge_engine.ErrorModelBResolved(); - error_model_a = em.a; - error_model_b = em.b; - } + // The stills scale and merge of --mode mx, on the frames as they were written. + auto sm = StillsScaleMerge(experiment, reflections, experiment.GetUnitCell(), dataset->stills_run_width, + reference_data, /*for_search=*/false, /*verdicts=*/false, nthreads, logger, + [](const std::string &) {}); + merged_reflections = std::move(sm.merged); + merged_statistics = std::move(sm.statistics); + reference_statistics = std::move(sm.reference_statistics); + resolution_fit_A = sm.resolution_fit_A; + error_model_isa = sm.isa; + error_model_isa_resolved = sm.isa_resolved; + error_model_a = sm.a; + error_model_b = sm.b; } logger.Info("Scale + merge completed in {:.2f} s ({} unique reflections)", diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 02364ce78..6fe078401 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -2527,18 +2527,22 @@ static Reflection MakeTestReflection(int i, int j) { .I = 1000.0f + f, .bkg = 10.0f + f, .sigma = 2.0f + 0.5f * f, + .dist_ewald = -1.5e-3f + 1e-4f * f, .prescaling_corr = 1.0f + 0.125f * f, .qe_corr = 1.0f / (1.0f + 0.01f * f), // always in (0, 1], as the correction itself is .flight_corr = 1.0f + 0.02f * f, // always >= 1, as the correction itself is .partiality = 0.5f + 0.01f * f, .zeta = 0.01f + 0.001f * f, - .image_scale_corr = 1.0f + 0.25f * f + .image_scale_corr = 1.0f + 0.25f * f, + .observed = j == 0, + .clipped = j == 1, + .overloaded = i == 3 }; } -// Every field of the round trip, against the reflection that was written. Deliberately NOT checked, -// because they are not part of it: dist_ewald and observed are prediction/integration scratch that -// is never written, and on_ice_ring is recomputed from the resolution by whoever scales (see the +// Every field of the round trip, against the reflection that was written, bit for bit: --mode scale merges +// what is read back, and reproduces the merge of the run that wrote it only if every field comes back as it +// was. Deliberately NOT checked: on_ice_ring, recomputed from the resolution by whoever scales (see the // ice-ring handling in Rugnux and in the --scale path). static void CheckReflectionRoundTrip(const Reflection &got, int i, int j) { INFO("image " << i << " reflection " << j); @@ -2546,22 +2550,26 @@ static void CheckReflectionRoundTrip(const Reflection &got, int i, int j) { CHECK(got.h == want.h); CHECK(got.k == want.k); CHECK(got.l == want.l); - CHECK(got.image_number == Catch::Approx(want.image_number)); - CHECK(got.delta_phi_deg == Catch::Approx(want.delta_phi_deg)); - CHECK(got.predicted_x == Catch::Approx(want.predicted_x)); - CHECK(got.predicted_y == Catch::Approx(want.predicted_y)); - CHECK(got.observed_x == Catch::Approx(want.observed_x)); - CHECK(got.observed_y == Catch::Approx(want.observed_y)); - CHECK(got.d == Catch::Approx(want.d)); - CHECK(got.I == Catch::Approx(want.I)); - CHECK(got.bkg == Catch::Approx(want.bkg)); - CHECK(got.sigma == Catch::Approx(want.sigma)); - CHECK(got.prescaling_corr == Catch::Approx(want.prescaling_corr)); // stored as 1/prescaling_corr, inverted again on read - CHECK(got.qe_corr == Catch::Approx(want.qe_corr)); // likewise stored as its reciprocal - CHECK(got.flight_corr == Catch::Approx(want.flight_corr)); // likewise - CHECK(got.partiality == Catch::Approx(want.partiality)); - CHECK(got.zeta == Catch::Approx(want.zeta)); - CHECK(got.image_scale_corr == Catch::Approx(want.image_scale_corr)); + CHECK(got.image_number == want.image_number); + CHECK(got.delta_phi_deg == want.delta_phi_deg); + CHECK(got.predicted_x == want.predicted_x); + CHECK(got.predicted_y == want.predicted_y); + CHECK(got.observed_x == want.observed_x); + CHECK(got.observed_y == want.observed_y); + CHECK(got.d == want.d); + CHECK(got.I == want.I); + CHECK(got.bkg == want.bkg); + CHECK(got.sigma == want.sigma); + CHECK(got.dist_ewald == want.dist_ewald); + CHECK(got.prescaling_corr == want.prescaling_corr); // stored as 1/prescaling_corr, inverted again on read + CHECK(got.qe_corr == want.qe_corr); // likewise stored as its reciprocal + CHECK(got.flight_corr == want.flight_corr); // likewise + CHECK(got.partiality == want.partiality); + CHECK(got.zeta == want.zeta); + CHECK(got.image_scale_corr == want.image_scale_corr); + CHECK(got.observed == want.observed); + CHECK(got.clipped == want.clipped); + CHECK(got.overloaded == want.overloaded); } TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") { @@ -2806,8 +2814,11 @@ TEST_CASE("JFJochReader_StillsRunWidth", "[HDF5][Full]") { REQUIRE_NOTHROW(writer.WriteHDF5(message)); EndMessage end_message; end_message.max_image_number = 1; - if (with_width) + if (with_width) { end_message.stills_run_width = w; + end_message.stills_given_cell = UnitCell{.a = 79.1f, .b = 79.1f, .c = 38.1f, + .alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f}; + } writer.WriteHDF5(end_message); writer.Finalize(); } @@ -2823,7 +2834,10 @@ TEST_CASE("JFJochReader_StillsRunWidth", "[HDF5][Full]") { CHECK(got->radial_scale == w.radial_scale); CHECK(got->table_dstar == w.table_dstar); CHECK(got->table_s == w.table_s); - } + REQUIRE(reader.GetDataset()->stills_given_cell.has_value()); + CHECK(reader.GetDataset()->stills_given_cell->c == 38.1f); + } else + CHECK_FALSE(reader.GetDataset()->stills_given_cell.has_value()); } remove("read_stills_run_width_master.h5"); remove("read_stills_run_width_data_000001.h5"); diff --git a/writer/HDF5DataFilePluginReflection.cpp b/writer/HDF5DataFilePluginReflection.cpp index 8ecfe7ee9..9e76f7e11 100644 --- a/writer/HDF5DataFilePluginReflection.cpp +++ b/writer/HDF5DataFilePluginReflection.cpp @@ -15,9 +15,13 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ return; std::vector h, k, l; - std::vector I, sigma, d, lp, qe, flight; + std::vector I, sigma, d; + // The three corrections are stored as their reciprocals, in double so that the reciprocal of the + // value read back is the float the reflection held. + std::vector lp, qe, flight; std::vector image, phi, pred_x, pred_y, obs_x, obs_y, bkg, var_bkg, partiality, zeta, - scale_factor; + scale_factor, dist_ewald; + std::vector observed, clipped, overloaded; h.reserve(msg.reflections.size()); k.reserve(msg.reflections.size()); @@ -39,6 +43,10 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ phi.reserve(msg.reflections.size()); zeta.reserve(msg.reflections.size()); scale_factor.reserve(msg.reflections.size()); + dist_ewald.reserve(msg.reflections.size()); + observed.reserve(msg.reflections.size()); + clipped.reserve(msg.reflections.size()); + overloaded.reserve(msg.reflections.size()); for (const auto &refl : msg.reflections) { image.emplace_back(refl.image_number); @@ -61,6 +69,10 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ phi.emplace_back(refl.delta_phi_deg); zeta.emplace_back(refl.zeta); scale_factor.emplace_back(refl.image_scale_corr); + dist_ewald.emplace_back(refl.dist_ewald); + observed.emplace_back(refl.observed); + clipped.emplace_back(refl.clipped); + overloaded.emplace_back(refl.overloaded); } std::string image_group_name = fmt::format("image_{:06d}", image_number); @@ -86,6 +98,10 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ image_group.SaveVector("partiality", partiality); image_group.SaveVector("zeta", zeta); image_group.SaveVector("image_scale_corr", scale_factor); + image_group.SaveVector("excitation_error", dist_ewald)->Units("Angstrom^-1"); + image_group.SaveVector("observed", observed); + image_group.SaveVector("clipped", clipped); + image_group.SaveVector("overloaded", overloaded); } void HDF5DataFilePluginReflection::WriteFinal(HDF5File &data_file) { diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 8836e4f96..0fcb295b4 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1208,6 +1208,11 @@ void NXmx::Finalize(const EndMessage &end) { SaveVector(*hdf5_file, "/entry/MX/stillsRockingWidthTable", w.table_s)->Units("Angstrom^-1"); } } + if (end.stills_given_cell) { + const auto &c = *end.stills_given_cell; + SaveVector(*hdf5_file, "/entry/MX/stillsGivenUnitCell", + std::vector{c.a, c.b, c.c, c.alpha, c.beta, c.gamma}); + } if (end.indexing_rate) { SaveScalar(*hdf5_file, "/entry/MX/imageIndexedMean", end.indexing_rate.value());