mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity
Files
fc6b91555f2d674e6b0c39e81154a66a1014ef09
Jungfraujoch/tests/WriteReflectionsTest.cpp
T
leonarski_fandClaude Opus 5.5 fc6b91555f Merged mmCIF: deposition keys, R_merge/R_pim/CC*, experiment metadata, scaled unmerged block 
The .cif now follows the PDBx/mmCIF dictionary's keys: _entry.id, the
entry_id of _cell/_symmetry/_reflns, pdbx_ordinal and pdbx_diffrn_id on
_reflns and _reflns_shell, diffrn_id on _diffrn_source/_diffrn_detector/
_diffrn_radiation, and _diffrn_radiation_wavelength.id. _software is a
loop naming Rugnux for 'data reduction' and for 'data scaling'.

The experiment as the input states it: _diffrn.ambient_temp from the
sample temperature, the source's general class (from NXmx
/entry/source/type, which the HDF5 reader now reads), facility and
beamline, the detector make/model in _diffrn_detector.type (where
.detector used to carry it) and its class in .detector where the name
settles it, and the collection date from /entry/start_time.

R_merge and R_pim are formed beside R_meas from the same per-reflection
sums (stills and rotation merges) and written with CC* overall and per
shell (_reflns.pdbx_Rmerge_I_obs/pdbx_Rpim_I_all/pdbx_CC_star and the
_reflns_shell counterparts), as report keys R_MERGE/R_PIM/CC_STAR and as
three new columns of the shell tables. R_pim following Weiss (2001),
CC* following Karplus & Diederichs (2012).

A second data block, data_sample_unmerged, carries the scaled unmerged
observations of a rotation merge as _diffrn_refln (measured index,
frame, scaled I and sigma, pdbx_scale_value) - the same records as the
.hkl, which the wwPDB accepts as scaled unmerged data. ScaledFull now
carries the scale and the frame. The unmerged MTZ stays unscaled.

Reflection data are unchanged: the merged _refln rows, the .mtz columns
and the .hkl are written as before.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
2026-10-09 17:50:11 +02:00

306 lines
14 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include <filesystem>
#include <fstream>
#include <string>
#include <cstring>
#include <gemmi/mtz.hpp>
#include <gemmi/read_cif.hpp>
#include <gemmi/numb.hpp>
#include "../common/DiffractionExperiment.h"
#include "../image_analysis/WriteReflections.h"
#include "../image_analysis/IntegrationOutcome.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include <cmath>
#include "SyntheticMergedReflections.h"
namespace {
// The wavelength CCP4's mtzlib substitutes when the data columns belong to dataset 0, which it
// takes for the reserved HKL_base: Cu K-alpha, and the wrong edge for anything reading f'/f''
// out of the file.
constexpr double CU_KALPHA_A = 1.54187;
constexpr UnitCell TETRAGONAL_CELL{47.0f, 47.0f, 63.0f, 90.0f, 90.0f, 90.0f};
DiffractionExperiment TestExperiment() {
DiffractionExperiment x;
x.IncidentEnergy_keV(12.7f); // ~0.976 A, nowhere near the Cu K-alpha default
x.SpaceGroupNumber(96); // P 43 21 2
x.SetUnitCell(TETRAGONAL_CELL);
return x;
}
}
TEST_CASE("Merged MTZ: the data dataset is id 1 and carries the wavelength", "[write_reflections][portable]") {
jfjoch_test::SyntheticMergeParams params;
params.true_space_group = "P 43 21 2";
params.twin_supergroup = "P 43 21 2";
params.d_min_A = 5.0;
const auto reflections = jfjoch_test::GenerateSyntheticMerged(params);
REQUIRE(!reflections.empty());
const auto experiment = TestExperiment();
const auto path = (std::filesystem::temp_directory_path() / "rugnux_merged_wavelength.mtz").string();
WriteMtzReflections(reflections, TETRAGONAL_CELL, experiment, path);
const gemmi::Mtz mtz = gemmi::read_mtz_file(path);
std::filesystem::remove(path);
// HKL_base at id 0, the data at id 1. A data dataset written at id 0 occupies the id MTZ
// reserves for the base, and mtzlib then reports CU_KALPHA_A instead of the real wavelength.
REQUIRE(mtz.datasets.size() == 2);
CHECK(mtz.datasets[0].id == 0);
CHECK(mtz.datasets[0].dataset_name == "HKL_base");
CHECK(mtz.datasets[1].id == 1);
CHECK(mtz.datasets[1].wavelength == Catch::Approx(experiment.GetWavelength_A()).epsilon(1e-5));
CHECK(mtz.datasets[1].wavelength != Catch::Approx(CU_KALPHA_A).epsilon(1e-3));
// The wavelength is read off the dataset the data columns belong to, so they have to be on the
// data dataset and not on the base.
for (const char *label : {"IMEAN", "SIGIMEAN", "F", "SIGF", "FreeR_flag"}) {
const gemmi::Mtz::Column *col = mtz.column_with_label(label);
REQUIRE(col != nullptr);
CHECK(col->dataset_id == 1);
}
// Cell and space group travel in the same header.
CHECK(mtz.spacegroup != nullptr);
CHECK(mtz.spacegroup->number == 96);
CHECK(mtz.cell.a == Catch::Approx(TETRAGONAL_CELL.a).epsilon(1e-5));
CHECK(mtz.cell.c == Catch::Approx(TETRAGONAL_CELL.c).epsilon(1e-5));
CHECK(mtz.cell.gamma == Catch::Approx(90.0).epsilon(1e-5));
CHECK(mtz.datasets[1].cell.a == Catch::Approx(TETRAGONAL_CELL.a).epsilon(1e-5));
CHECK(mtz.nreflections == static_cast<int>(reflections.size()));
}
TEST_CASE("Unmerged MTZ: LP is Lorentz-polarization, QE the sensor efficiency, FLIGHT the flight path",
"[write_reflections][portable]") {
// The whole point of the split: LP must mean what XDS and DIALS mean by it, and the raw count
// sum must still be recoverable from the file alone, as I / LP * QE * FLIGHT.
auto experiment = TestExperiment();
experiment.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}));
IntegrationOutcome outcome;
const float raw[3] = {1000.0f, 250.0f, 40.0f};
const float lp[3] = {1.75f, 2.50f, 0.90f}; // Lorentz x polarization, and nothing else
const float qe[3] = {0.9375f, 0.8125f, 1.0f}; // 1.0 = the sensor said nothing to correct
const float fl[3] = {1.0125f, 1.0400f, 1.0f}; // >= 1: an oblique reflection crossed more air
for (int i = 0; i < 3; ++i) {
Reflection r{};
r.h = 4 + i; r.k = 2; r.l = 6;
r.image_number = static_cast<float>(i);
r.d = 5.0f + i;
r.I = raw[i]; // the writer is what applies the factor
r.sigma = std::sqrt(raw[i]);
r.prescaling_corr = lp[i];
r.qe_corr = qe[i];
r.flight_corr = fl[i];
r.partiality = 1.0f;
r.predicted_x = 100.0f + i; r.predicted_y = 200.0f + i;
r.observed_x = NAN; r.observed_y = NAN;
outcome.reflections.push_back(r);
}
const auto path = (std::filesystem::temp_directory_path() / "rugnux_unmerged_qe.mtz").string();
WriteUnmergedMtzReflections({outcome}, TETRAGONAL_CELL, experiment, false, path);
const gemmi::Mtz mtz = gemmi::read_mtz_file(path);
std::filesystem::remove(path);
const gemmi::Mtz::Column *c_I = mtz.column_with_label("I");
const gemmi::Mtz::Column *c_lp = mtz.column_with_label("LP");
const gemmi::Mtz::Column *c_qe = mtz.column_with_label("QE");
const gemmi::Mtz::Column *c_fl = mtz.column_with_label("FLIGHT");
REQUIRE(c_I != nullptr);
REQUIRE(c_lp != nullptr);
REQUIRE(c_qe != nullptr); // DIALS writes this column even when there is nothing in it
REQUIRE(c_fl != nullptr);
REQUIRE(mtz.nreflections == 3);
for (int i = 0; i < 3; ++i) {
const float I = mtz.data[i * mtz.columns.size() + c_I->idx];
const float LP = mtz.data[i * mtz.columns.size() + c_lp->idx];
const float QE = mtz.data[i * mtz.columns.size() + c_qe->idx];
const float FL = mtz.data[i * mtz.columns.size() + c_fl->idx];
INFO("row " << i);
// LP holds Lorentz x polarization alone: the sensor term was never inside it.
CHECK(LP == Catch::Approx(lp[i]).epsilon(1e-5));
// QE is a divisor normalised to 1 at normal incidence, so it never drops below 1.
CHECK(QE == Catch::Approx(1.0f / qe[i]).epsilon(1e-5));
CHECK(QE >= 1.0f);
// FLIGHT is a divisor in the same convention, and runs the other way: the sensor favours an
// oblique reflection, the medium attenuates it, so this one never rises above 1.
CHECK(FL == Catch::Approx(1.0f / fl[i]).epsilon(1e-5));
CHECK(FL <= 1.0f);
// ... and the three together put the raw counts back.
CHECK(I / LP * QE * FL == Catch::Approx(raw[i]).epsilon(1e-4));
// The intensity itself is the fully corrected value - all three applied.
CHECK(I == Catch::Approx(raw[i] * lp[i] * qe[i] * fl[i]).epsilon(1e-5));
}
}
TEST_CASE("Unmerged MTZ: built on several workers, the same rows as on one", "[write_reflections][portable]") {
// A sweep large enough that the sort, the event sums and the row fill are all split into
// pieces, with events of several parts, events cut short by the sweep's ends, and parts that
// share (h,k,l) and image_number - the tie the sort has to break the same way every time.
auto experiment = TestExperiment();
experiment.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}));
uint32_t seed = 12345;
const auto next = [&seed] { seed = seed * 1664525u + 1013904223u; return (seed >> 8) / 16777216.0f; };
std::vector<IntegrationOutcome> outcomes(200);
for (int f = 0; f < 200; ++f) {
for (int q = 0; q < 8000; ++q) {
if ((f + q) % 7 >= 3)
continue; // each reflection is seen on three consecutive images out of seven
Reflection r{};
r.h = q % 20 - 10; r.k = (q / 20) % 20 - 10; r.l = q / 400 + 1;
r.image_number = static_cast<float>(f);
r.I = 1000.0f * next();
r.sigma = 1.0f + 10.0f * next();
r.bkg = next();
r.var_bkg = next();
r.prescaling_corr = 1.0f + next();
r.partiality = 0.2f + 0.2f * next();
r.predicted_x = 1000.0f * next(); r.predicted_y = 1000.0f * next();
r.observed_x = NAN; r.observed_y = NAN;
outcomes[f].reflections.push_back(r);
if (q % 97 == 0) {
r.I = 1000.0f * next();
outcomes[f].reflections.push_back(r);
}
}
}
for (bool sum_partials : {true, false}) {
INFO("sum_partials " << sum_partials);
const gemmi::Mtz serial = UnmergedMtz(outcomes, TETRAGONAL_CELL, experiment, sum_partials, 1);
const gemmi::Mtz parallel = UnmergedMtz(outcomes, TETRAGONAL_CELL, experiment, sum_partials, 8);
REQUIRE(serial.nreflections > 100000);
REQUIRE(parallel.nreflections == serial.nreflections);
CHECK(std::memcmp(parallel.data.data(), serial.data.data(), serial.data.size() * sizeof(float)) == 0);
REQUIRE(parallel.batches.size() == serial.batches.size());
CHECK(parallel.batches.front().number == serial.batches.front().number);
CHECK(parallel.batches.back().number == serial.batches.back().number);
}
}
TEST_CASE("SHELX .hkl of unmerged fulls: fixed 3I4,2F8.2, scaled to fit, 0 0 0 terminator",
"[write_reflections][portable]") {
// Two fulls of the same reflection and a Friedel mate: written as measured, nothing averaged. The
// largest |I| or sigma is put at 9999.00; an unusable full (sigma not positive) is left out.
const std::vector<ScaledFull> fulls{
{1, 2, 3, 200000.0f, 2000.0f},
{-1, -2, -3, 100000.0f, 1200.0f},
{1, 2, 3, -50.0f, 20.0f},
{4, 5, 6, 10.0f, 0.0f},
};
const auto path = (std::filesystem::temp_directory_path() / "jfjoch_scaled_fulls_test.hkl").string();
WriteShelxHklReflections(fulls, path, 2);
std::vector<std::string> lines;
{
// Closed before the remove: Windows refuses to delete a file that is still open.
std::ifstream in(path);
for (std::string line; std::getline(in, line);)
lines.push_back(line);
}
std::filesystem::remove(path);
REQUIRE(lines.size() == 4);
for (const auto &line : lines)
CHECK(line.size() == 28);
CHECK(lines[0] == " 1 2 3 9999.00 99.99");
CHECK(lines[1] == " -1 -2 -3 4999.50 59.99");
CHECK(lines[2] == " 1 2 3 -2.50 1.00");
CHECK(lines[3] == " 0 0 0 0.00 0.00");
}
TEST_CASE("CC* from CC1/2", "[write_reflections][portable]") {
CHECK(CCStar(1.0) == Catch::Approx(1.0));
CHECK(CCStar(0.5) == Catch::Approx(std::sqrt(2.0 / 3.0)));
CHECK(std::isnan(CCStar(0.0)));
CHECK(std::isnan(CCStar(-0.1)));
CHECK(std::isnan(CCStar(NAN)));
}
TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and the scaled unmerged block",
"[write_reflections][portable]") {
jfjoch_test::SyntheticMergeParams params;
params.true_space_group = "P 43 21 2";
params.twin_supergroup = "P 43 21 2";
params.d_min_A = 5.0;
const auto reflections = jfjoch_test::GenerateSyntheticMerged(params);
REQUIRE(!reflections.empty());
auto experiment = TestExperiment();
experiment.SampleTemperature_K(100.0f);
MergeStatistics statistics;
MergeStatisticsShell shell;
shell.d_min = 5.0f; shell.d_max = 30.0f;
shell.total_observations = 400; shell.unique_reflections = 100; shell.possible_unique_reflections = 110;
shell.mean_i_over_sigma = 10.0;
shell.cc_half = 0.98;
shell.r_merge = 0.05; shell.r_meas = 0.06; shell.r_pim = 0.03;
statistics.shells = {shell};
statistics.overall = shell;
// Three fulls, one without a usable sigma - that one is left out, as from the .hkl.
const std::vector<ScaledFull> fulls{
{1, 2, 3, 200.0f, 20.0f, 1.25f, 0.4f},
{-1, -2, -3, 100.0f, 12.0f, 0.8f, 9.6f},
{4, 5, 6, 10.0f, 0.0f, 1.0f, 3.0f},
};
OutputProvenance provenance;
provenance.collection_date = "2026-07-13T13:48:52.413Z";
const auto path = (std::filesystem::temp_directory_path() / "rugnux_merged_test.cif").string();
WriteMmcifReflections(reflections, TETRAGONAL_CELL, experiment, statistics, ErrorModelReport{},
TwinningAnalysisResult{}, fulls, provenance, path, 2);
gemmi::cif::Document doc = gemmi::read_cif_gz(path);
std::filesystem::remove(path);
REQUIRE(doc.blocks.size() == 2);
gemmi::cif::Block &merged = doc.blocks[0];
CHECK(merged.name == "sample");
CHECK(*merged.find_value("_reflns.pdbx_ordinal") == "1");
CHECK(*merged.find_value("_reflns.pdbx_diffrn_id") == "1");
CHECK(*merged.find_value("_reflns.pdbx_Rmerge_I_obs") == "0.0500");
CHECK(*merged.find_value("_reflns.pdbx_Rpim_I_all") == "0.0300");
CHECK(gemmi::cif::as_number(*merged.find_value("_reflns.pdbx_CC_star"))
== Catch::Approx(std::sqrt(2 * 0.98 / 1.98)).margin(1e-4));
CHECK(*merged.find_value("_diffrn.ambient_temp") == "100.0");
CHECK(*merged.find_value("_diffrn_detector.pdbx_collection_date") == "2026-07-13");
gemmi::cif::Table software = merged.find("_software.", {"pdbx_ordinal", "classification"});
REQUIRE(software.length() == 2);
CHECK(gemmi::cif::as_string(software[0][1]) == "data reduction");
CHECK(gemmi::cif::as_string(software[1][1]) == "data scaling");
gemmi::cif::Table shells = merged.find("_reflns_shell.", {"pdbx_ordinal", "Rmerge_I_obs",
"pdbx_Rpim_I_all", "pdbx_CC_star"});
REQUIRE(shells.length() == 1);
CHECK(shells[0][1] == "0.0500");
CHECK(shells[0][2] == "0.0300");
gemmi::cif::Block &unmerged = doc.blocks[1];
CHECK(unmerged.name == "sample_unmerged");
gemmi::cif::Table obs = unmerged.find("_diffrn_refln.", {"id", "pdbx_image_id", "index_h",
"intensity_net", "intensity_sigma",
"pdbx_scale_value"});
REQUIRE(obs.length() == 2);
CHECK(obs[0][0] == "1");
CHECK(obs[0][1] == "1"); // peaks on frame 0.4 -> image 1
CHECK(obs[1][1] == "11"); // 9.6 -> image 11
CHECK(obs[1][2] == "-1");
CHECK(gemmi::cif::as_number(obs[0][3]) == Catch::Approx(200.0));
CHECK(gemmi::cif::as_number(obs[1][5]) == Catch::Approx(0.8));
}
Reference in New Issue View Git Blame Copy Permalink