Merge stills-bugs into stills-combined: --wavelength print, _process.h5 without image data, shared StillsScaleMerge for mx and --mode scale
# Conflicts: # rugnux/RugnuxScaleMerge.cpp
This commit is contained in:
leonarski_f committed
commit
199f5919ee
15 files changed
+514
-351
No files matched your search
@@ -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<std::array<int32_t, 9>> 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<StillsRunWidth> stills_run_width;
|
||||
// Stills: the cell given on the command line, whose population of frames the merge keeps.
|
||||
std::optional<UnitCell> 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
|
||||
|
||||
+5
-2
@@ -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
|
||||
|
||||
+10
-7
@@ -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
|
||||
|
||||
+103
-76
@@ -331,18 +331,24 @@ bool ReadReflectionsFromGroup(HDF5Object &file,
|
||||
auto bkg = file.ReadOptVector<float>(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<float>(image_group_name + "/background_variance");
|
||||
auto lp = file.ReadOptVector<float>(image_group_name + "/lp");
|
||||
auto lp = file.ReadOptVector<double>(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<float>(image_group_name + "/qe");
|
||||
auto qe = file.ReadOptVector<double>(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<float>(image_group_name + "/flight");
|
||||
auto flight = file.ReadOptVector<double>(image_group_name + "/flight");
|
||||
auto partiality = file.ReadOptVector<float>(image_group_name + "/partiality");
|
||||
auto phi = file.ReadOptVector<float>(image_group_name + "/delta_phi");
|
||||
auto zeta = file.ReadOptVector<float>(image_group_name + "/zeta");
|
||||
auto image_scale_corr = file.ReadOptVector<float>(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<float>(image_group_name + "/excitation_error");
|
||||
auto observed = file.ReadOptVector<uint8_t>(image_group_name + "/observed");
|
||||
auto clipped = file.ReadOptVector<uint8_t>(image_group_name + "/clipped");
|
||||
auto overloaded = file.ReadOptVector<uint8_t>(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<float>(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<float>(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<float>(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<float>(
|
||||
"/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image");
|
||||
else
|
||||
dataset->efficiency = std::vector<float>(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<int64_t>("/entry/roi/" + s + "/max"));
|
||||
dataset->roi_sum.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/sum"));
|
||||
dataset->roi_sum_sq.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/sum_sq"));
|
||||
dataset->roi_npixel.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/npixel"));
|
||||
dataset->roi_x.emplace_back(master_file->ReadVector<float>("/entry/roi/" + s + "/x"));
|
||||
dataset->roi_y.emplace_back(master_file->ReadVector<float>("/entry/roi/" + s + "/y"));
|
||||
}
|
||||
|
||||
if (master_file->Exists("/entry/MX")) {
|
||||
if (master_file->Exists("/entry/MX/peakCountUnfiltered"))
|
||||
dataset->spot_count = master_file->ReadOptVector<float>("/entry/MX/peakCountUnfiltered");
|
||||
else
|
||||
dataset->spot_count = master_file->ReadOptVector<float>("/entry/MX/nPeaks");
|
||||
|
||||
dataset->spot_count_low_res = master_file->ReadOptVector<float>("/entry/MX/peakCountLowRes");
|
||||
dataset->spot_count_indexed = master_file->ReadOptVector<float>("/entry/MX/peakCountIndexed");
|
||||
dataset->spot_count_ice_rings = master_file->ReadOptVector<float>("/entry/MX/peakCountIceRingRes");
|
||||
dataset->spot_count_ice_control = master_file->ReadOptVector<float>("/entry/MX/peakCountIceRingControl");
|
||||
|
||||
dataset->indexing_result = master_file->ReadOptVector<float>("/entry/MX/imageIndexed");
|
||||
dataset->bkg_estimate = master_file->ReadOptVector<float>("/entry/MX/bkgEstimate");
|
||||
dataset->spindle_blind_fraction = master_file->ReadOptVector<float>("/entry/MX/spindleBlindFraction");
|
||||
dataset->ice_ring_score = master_file->ReadOptVector<float>("/entry/MX/iceRingScore");
|
||||
dataset->resolution_estimate = master_file->ReadOptVector<float>("/entry/MX/resolutionEstimate");
|
||||
dataset->profile_radius = master_file->ReadOptVector<float>("/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<float>("/entry/MX/indexedLatticeCount");
|
||||
if (dataset->indexing_lattice_count.empty())
|
||||
dataset->indexing_lattice_count = master_file->ReadOptVector<float>("/entry/MX/indexingLatticeCount");
|
||||
dataset->mosaicity_deg = master_file->ReadOptVector<float>("/entry/MX/mosaicity");
|
||||
dataset->b_factor = master_file->ReadOptVector<float>("/entry/MX/bFactor");
|
||||
dataset->image_scale_factor = master_file->ReadOptVector<float>("/entry/MX/imageScaleFactor");
|
||||
dataset->image_scale_cc = master_file->ReadOptVector<float>("/entry/MX/imageScaleCC");
|
||||
dataset->integrated_reflections = master_file->ReadOptVector<float>("/entry/MX/integratedReflections");
|
||||
dataset->sweep_quality = master_file->ReadOptVector<uint8_t>("/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<std::string>("/entry/MX/sweepQualityReasons", i)
|
||||
.value_or(""));
|
||||
}
|
||||
dataset->frame_disposition = master_file->ReadOptVector<uint8_t>("/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<std::string>("/entry/MX/frameDispositionCodes", i)
|
||||
.value_or(""));
|
||||
}
|
||||
}
|
||||
if (master_file->Exists("/entry/image"))
|
||||
dataset->max_value = master_file->ReadOptVector<int64_t>("/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<float>(
|
||||
"/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image");
|
||||
else
|
||||
dataset->efficiency = std::vector<float>(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<int64_t>("/entry/roi/" + s + "/max"));
|
||||
dataset->roi_sum.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/sum"));
|
||||
dataset->roi_sum_sq.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/sum_sq"));
|
||||
dataset->roi_npixel.emplace_back(master_file->ReadVector<int64_t>("/entry/roi/" + s + "/npixel"));
|
||||
dataset->roi_x.emplace_back(master_file->ReadVector<float>("/entry/roi/" + s + "/x"));
|
||||
dataset->roi_y.emplace_back(master_file->ReadVector<float>("/entry/roi/" + s + "/y"));
|
||||
}
|
||||
|
||||
if (master_file->Exists("/entry/MX")) {
|
||||
if (master_file->Exists("/entry/MX/peakCountUnfiltered"))
|
||||
dataset->spot_count = master_file->ReadOptVector<float>("/entry/MX/peakCountUnfiltered");
|
||||
else
|
||||
dataset->spot_count = master_file->ReadOptVector<float>("/entry/MX/nPeaks");
|
||||
|
||||
dataset->spot_count_low_res = master_file->ReadOptVector<float>("/entry/MX/peakCountLowRes");
|
||||
dataset->spot_count_indexed = master_file->ReadOptVector<float>("/entry/MX/peakCountIndexed");
|
||||
dataset->spot_count_ice_rings = master_file->ReadOptVector<float>("/entry/MX/peakCountIceRingRes");
|
||||
dataset->spot_count_ice_control = master_file->ReadOptVector<float>("/entry/MX/peakCountIceRingControl");
|
||||
|
||||
dataset->indexing_result = master_file->ReadOptVector<float>("/entry/MX/imageIndexed");
|
||||
dataset->bkg_estimate = master_file->ReadOptVector<float>("/entry/MX/bkgEstimate");
|
||||
dataset->spindle_blind_fraction = master_file->ReadOptVector<float>("/entry/MX/spindleBlindFraction");
|
||||
dataset->ice_ring_score = master_file->ReadOptVector<float>("/entry/MX/iceRingScore");
|
||||
dataset->resolution_estimate = master_file->ReadOptVector<float>("/entry/MX/resolutionEstimate");
|
||||
dataset->profile_radius = master_file->ReadOptVector<float>("/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<float>("/entry/MX/indexedLatticeCount");
|
||||
if (dataset->indexing_lattice_count.empty())
|
||||
dataset->indexing_lattice_count = master_file->ReadOptVector<float>("/entry/MX/indexingLatticeCount");
|
||||
dataset->mosaicity_deg = master_file->ReadOptVector<float>("/entry/MX/mosaicity");
|
||||
dataset->b_factor = master_file->ReadOptVector<float>("/entry/MX/bFactor");
|
||||
dataset->image_scale_factor = master_file->ReadOptVector<float>("/entry/MX/imageScaleFactor");
|
||||
dataset->image_scale_cc = master_file->ReadOptVector<float>("/entry/MX/imageScaleCC");
|
||||
dataset->integrated_reflections = master_file->ReadOptVector<float>("/entry/MX/integratedReflections");
|
||||
dataset->sweep_quality = master_file->ReadOptVector<uint8_t>("/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<std::string>("/entry/MX/sweepQualityReasons", i)
|
||||
.value_or(""));
|
||||
}
|
||||
dataset->frame_disposition = master_file->ReadOptVector<uint8_t>("/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<std::string>("/entry/MX/frameDispositionCodes", i)
|
||||
.value_or(""));
|
||||
}
|
||||
}
|
||||
if (master_file->Exists("/entry/image"))
|
||||
dataset->max_value = master_file->ReadOptVector<int64_t>("/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<float>("/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<IntegrationOutcome> HDF5MetadataSource::ReadReflections(size_t start
|
||||
// ── per-image mosaicity ───────────────────────────────────────────────
|
||||
if (meta_file->Exists("/entry/MX/mosaicity")) {
|
||||
try {
|
||||
outcome.mosaicity_deg =
|
||||
meta_file->ReadElement<float>("/entry/MX/mosaicity", meta_image_id);
|
||||
// NaN is how the file says the image had none.
|
||||
const auto mos = meta_file->ReadElement<float>("/entry/MX/mosaicity", meta_image_id);
|
||||
if (mos && std::isfinite(*mos))
|
||||
outcome.mosaicity_deg = mos;
|
||||
} catch (...) {
|
||||
}
|
||||
}
|
||||
|
||||
@@ -15,6 +15,7 @@ void JFJochHDF5Reader::ReadFile(const std::string &filename) {
|
||||
auto metadata = std::make_shared<HDF5MetadataSource>();
|
||||
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<JFJochReaderDataset> &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,
|
||||
|
||||
@@ -31,6 +31,9 @@ class JFJochHDF5Reader : public JFJochReader {
|
||||
std::shared_ptr<HDF5MetadataSource> 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<JFJochReaderDataset> &dataset,
|
||||
DataMessage &message,
|
||||
|
||||
@@ -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<UnitCell> 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<StillsRunWidth> stills_run_width;
|
||||
// Stills: the cell the run was given, whose population of frames it merged (/entry/MX/stillsGivenUnitCell).
|
||||
std::optional<UnitCell> stills_given_cell;
|
||||
|
||||
std::vector<float> az_int_bin_to_q;
|
||||
std::vector<float> az_int_bin_to_phi;
|
||||
|
||||
+6
-8
@@ -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<JFJochHDF5Reader *>(&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");
|
||||
|
||||
+31
-1
@@ -356,6 +356,34 @@ bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const Vali
|
||||
std::pair<double, double> 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<IntegrationOutcome> &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<MergedReflection> merged;
|
||||
MergeStatistics statistics;
|
||||
std::optional<MergeStatistics> reference_statistics; // --report-resolution
|
||||
std::optional<double> d_min; // the cut applied to merged
|
||||
std::optional<double> resolution_fit_A;
|
||||
std::optional<StillsRunWidth> used_width; // the curve the post-refinement used
|
||||
std::vector<std::vector<ObservationVerdict>> 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<IntegrationOutcome> &outcomes,
|
||||
const std::optional<UnitCell> &reference_cell,
|
||||
const std::optional<StillsRunWidth> &run_width,
|
||||
const std::vector<MergedReflection> &reference_data, bool for_search,
|
||||
bool verdicts, int nthreads, Logger &logger,
|
||||
const std::function<void(const std::string &)> &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<StillsPartialityRefine::RunWidth> 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<std::vector<ObservationVerdict>> observation_verdicts_;
|
||||
std::optional<double> observation_d_min_;
|
||||
std::optional<StillsRunWidth> 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).
|
||||
|
||||
+191
-148
@@ -599,6 +599,173 @@ static void WarnIfMetricHostsMoreSymmetry(const std::optional<UnitCell> &cell,
|
||||
metric_system, metric_rotations, sg.xhm(), group_rotations);
|
||||
}
|
||||
|
||||
size_t LeaveOutOtherCellForms(std::vector<IntegrationOutcome> &outcomes, const UnitCell &given_cell, Logger &logger) {
|
||||
const double given = CrystalLattice(given_cell).CalcVolume();
|
||||
std::vector<double> 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<double>::infinity();
|
||||
high = std::numeric_limits<double>::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<IntegrationOutcome> &outcomes,
|
||||
const std::optional<UnitCell> &reference_cell,
|
||||
const std::optional<StillsRunWidth> &run_width,
|
||||
const std::vector<MergedReflection> &reference_data, bool for_search,
|
||||
bool verdicts, int nthreads, Logger &logger,
|
||||
const std::function<void(const std::string &)> &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<float> cc, g; std::optional<int64_t> cc_n; std::vector<float> corr; };
|
||||
std::vector<FrameScale> 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<double> 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<double>::infinity();
|
||||
high = std::numeric_limits<double>::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<IntegrationOutcome> &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<double> 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()));
|
||||
|
||||
+51
-82
@@ -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<float>(r.h) + o.latt.Bstar() * static_cast<float>(r.k)
|
||||
+ o.latt.Cstar() * static_cast<float>(r.l) + S0;
|
||||
r.dist_ewald = static_cast<float>(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<double> 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)",
|
||||
|
||||
+36
-22
@@ -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");
|
||||
|
||||
@@ -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<std::unique_ptr<JFJochReader>> inputs;
|
||||
auto input = std::make_unique<JFJochHDF5Reader>();
|
||||
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);
|
||||
|
||||
|
||||
@@ -15,9 +15,13 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_
|
||||
return;
|
||||
|
||||
std::vector<int32_t> h, k, l;
|
||||
std::vector<float> I, sigma, d, lp, qe, flight;
|
||||
std::vector<float> 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<double> lp, qe, flight;
|
||||
std::vector<float> image, phi, pred_x, pred_y, obs_x, obs_y, bkg, var_bkg, partiality, zeta,
|
||||
scale_factor;
|
||||
scale_factor, dist_ewald;
|
||||
std::vector<uint8_t> 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) {
|
||||
|
||||
@@ -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<float>{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());
|
||||
|
||||
Reference in new issue
Block a user