mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions 2 Packages Projects Releases 72 Wiki Activity

_process.h5: record the change of basis, and the space group

Browse Source 
A _process.h5 was internally inconsistent whenever the space group re-seated
the lattice. The per-image reflections and lattices go to file as each image
is processed, but the conventional setting is only chosen afterwards, so the
file kept pre-reindex indices beside a post-reindex cell. Measured on an
affected dataset, max|d_file - d(hkl, stored cell)| was 18.6 A. --mode scale
then compared each frame's lattice against the stored cell, found none of
1800 acceptable, rejected every observation and died in the merge with
"resolution calculation failed".

/entry/MX/reindexMatrix now carries M with hkl_cell = M . hkl_written, and
the reader applies it, so everything it hands out is in the setting of
/entry/sample/unit_cell. Absent means the identity, so a file written before
this reads exactly as before. On the affected dataset M comes out
[[1,1,0],[0,1,1],[1,0,1]], det 2 - the primitive-to-body-centred basis its
volume ratio implied - and the same measure falls from 18.6 A to 1.8e-5.

Writing the reflections in the final setting instead was rejected: the
per-image writer is shared with the broker, which streams and cannot buffer a
run; and h,k,l, predicted_x/y and the per-image lattice are one consistent
statement about one image, which retro-editing the indices would silently
break.

Two things turned up while fixing it. There are three re-seat sites, not one
- the space-group search's own centred-lattice test re-seats too, and logged
nothing - so the matrix composes over all of them. And the space group itself
was never written: it was set only on the arm that searches, while a two-pass
rotation run reuses pass 1's group and takes the other arm, so the canonical
file carried a cell but no group and --mode scale merged in P1.

--mode scale now reproduces --mode mx on the affected dataset: same space
group and cell, 39329 unique reflections both, 100% of reflections common,
CC 0.99985, sum|dI|/sum|I| = 0.0081. The residual is three reflections in
11.09 M crossing an ice-band edge, because mx carries the integrator's d and
scale recomputes it from the cell. An unaffected dataset is byte-identical in
.mtz, .hkl and .cif.

An older affected file still cannot merge - M is not recoverable from it -
but now says so in 51 s, naming both cells and the -S/-C override to use,
instead of failing inside the merge.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
This commit is contained in:
leonarski_f
2026-08-25 21:21:55 +02:00
co-authored by Claude Opus 5
parent 8fa1ab7b05
commit db908c6ec1
10 changed files with 266 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/JFJochReaderTest.cpp
+111
View File
@@ -2479,6 +2479,117 @@ TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") {
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// The per-image reflections and lattices are written in the setting the images were indexed in, but
// the space group - and with it the conventional setting the cell beside them is in - is only settled
// after the merge, so the two can differ by an integral change of basis. /entry/MX/reindexMatrix
// carries it, and the reader applies it, so what comes out is in the cell's setting. A file without
// the dataset (every file written before it existed) is read as the identity, which is what the
// round-trip tests above check.
TEST_CASE("JFJochReader_ReadReflections_Reindex", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("read_reflections_reindex")
.ImagesPerTrigger(2)
.ImagesPerFile(1)
.OverwriteExistingFiles(true)
.BitDepthImage(16)
.PixelSigned(true)
.SetFileWriterFormat(FileWriterFormat::NXmxVDS)
.IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
.Compression(CompressionAlgorithm::NO_COMPRESSION);
// hkl_cell = M . hkl_written, det 2 - the size of step a primitive-to-centred re-seat takes.
const std::array<int32_t, 9> M = {1, 1, 0,
0, 1, 1,
1, 0, 1};
std::vector<int16_t> image(x.GetPixelsNum(), 0);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter writer(start_message);
ScanResultGenerator scan_result(x);
for (int i = 0; i < x.GetImageNum(); i++) {
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
if (i == 1) {
message.integrated_reflections = 2;
message.reflections = {MakeTestReflection(i, 0), MakeTestReflection(i, 1)};
message.indexing_result = true;
message.indexing_lattice = CrystalLattice({100,0,0}, {0,50,0}, {0,0,30});
}
REQUIRE_NOTHROW(writer.WriteHDF5(message));
scan_result.Add(message);
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
end_message.reindex_matrix = M;
scan_result.FillEndMessage(end_message);
writer.WriteHDF5(end_message);
writer.Finalize();
}
// hkl and the lattice come back in the cell's setting; every other field is untouched.
const auto check = [&](const Reflection &got, int j) {
const Reflection want = MakeTestReflection(1, j);
INFO("reflection " << j);
CHECK(got.h == M[0] * want.h + M[1] * want.k + M[2] * want.l);
CHECK(got.k == M[3] * want.h + M[4] * want.k + M[5] * want.l);
CHECK(got.l == M[6] * want.h + M[7] * want.k + M[8] * want.l);
CHECK(got.I == Catch::Approx(want.I));
CHECK(got.d == Catch::Approx(want.d));
CHECK(got.image_number == Catch::Approx(want.image_number));
};
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("read_reflections_reindex_master.h5"));
REQUIRE(reader.GetDataset()->reindex_matrix.has_value());
CHECK(reader.GetDataset()->reindex_matrix.value() == M);
auto reflections = reader.ReadReflections();
REQUIRE(reflections.size() == 2);
REQUIRE(reflections[1].reflections.size() == 2);
check(reflections[1].reflections[0], 0);
check(reflections[1].reflections[1], 1);
// latt = M . latt_written, row by row: (100,0,0)+(0,50,0), (0,50,0)+(0,0,30), (100,0,0)+(0,0,30).
CHECK(reflections[1].latt.Vec0().x == Catch::Approx(100.0f));
CHECK(reflections[1].latt.Vec0().y == Catch::Approx(50.0f));
CHECK(reflections[1].latt.Vec1().y == Catch::Approx(50.0f));
CHECK(reflections[1].latt.Vec1().z == Catch::Approx(30.0f));
CHECK(reflections[1].latt.Vec2().x == Catch::Approx(100.0f));
CHECK(reflections[1].latt.Vec2().z == Catch::Approx(30.0f));
CHECK(reflections[1].latt.CalcVolume() == Catch::Approx(2.0 * 100 * 50 * 30));
// The per-image message path (the viewer's) is re-seated the same way.
auto reader_image = reader.LoadImage(1);
REQUIRE(reader_image);
REQUIRE(reader_image->ImageData().reflections.size() == 2);
check(reader_image->ImageData().reflections[0], 0);
REQUIRE(reader_image->ImageData().indexing_lattice);
CHECK(reader_image->ImageData().indexing_lattice->CalcVolume()
== Catch::Approx(2.0 * 100 * 50 * 30));
}
remove("read_reflections_reindex_master.h5");
remove("read_reflections_reindex_data_000001.h5");
remove("read_reflections_reindex_data_000002.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
static std::vector<SpotToSave> MakeTestSpots(int i) {
return {
SpotToSave{