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 e6908eab6..cca40f5f7 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 @@ -397,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 0d6fa3b6d..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); } @@ -538,6 +548,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 +633,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 +810,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")) { @@ -973,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")); @@ -1676,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/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/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 10fd7e204..fbc18ffb2 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1035,11 +1035,11 @@ 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. + // 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. @@ -1063,8 +1063,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b ? -(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; - } else { - start_message.write_images = true; // no file to link to; the process file carries them } // The full-detector mask is 4 B/px and the ROI map 2 B/px - tens of megabytes each - and nothing @@ -1313,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 a6b1c4870..4d870b550 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -599,6 +599,173 @@ 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 (!for_search) + logger.Info("Merge: {} merged observations have a partiality below {}", merge_engine.FarTailCount(), + MergeOnTheFly::FAR_TAIL_PARTIALITY); + 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 +862,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,106 +1852,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 (!for_search && write_output) - logger.Info("Merge: {} merged observations have a partiality below {}", merge_engine.FarTailCount(), - MergeOnTheFly::FAR_TAIL_PARTIALITY); - 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 d551391bf..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" @@ -2163,15 +2164,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 +2189,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"))); } @@ -2433,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 @@ -2509,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() @@ -2542,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/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 7114e22f3..fa971397c 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_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_noimg_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_noimg_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); + + { + HDF5ReadOnlyFile file("process_noimg_out_process.h5"); + CHECK_FALSE(file.Exists("/entry/data/data")); + } + { + JFJochHDF5Reader out; + REQUIRE_NOTHROW(out.ReadFile("process_noimg_out_process.h5")); + CHECK(out.GetNumberOfImages() == 4); + CHECK(out.LoadImage(0) == nullptr); + CHECK_NOTHROW(out.ReadReflections(0, 3)); + } + + reader.Close(); + 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); +} + TEST_CASE("Rugnux_NoOutput", "[HDF5][Full]") { WriteTestDataset("process_noout_in", 6); 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());