mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 3 Packages Projects Releases 80 Wiki Activity
Files
eeeb7248906289931a80070dcdd89878b2aa0509
Jungfraujoch/tests/RugnuxTest.cpp
T
leonarski_fandClaude Opus 5.5 eeeb724890 rugnux: decide the two pass decisions on evidence net of its own noise 
Two decisions in the rotation two-pass flipped whenever the integration
changed, because each compared a raw number with no account of its noise.

The quality guard (RefinedPassIsWorse) counted reflections at I/sigma >= 2 in
each pass's search merge. Across two lattices that count is not like for like:
a header pass on a 7x supercell merges seven times the reflections, six in
seven empty, and their noise reached the strong count - it out-counted the
refined pass on the crystal's own lattice and the run went back to the
supercell, while that merge held more reflections at I/sigma <= -2 than at
>= 2. The guard now reads the signal as the strong count less the count at
<= -2 (the merge's own noise tail), and its dilution clause scales the
reflection count of the smaller lattice by the integer index of one lattice in
the other, so the same lattice (index 1) compares exactly as before. Measured
on 37 sets: no guard decision changes except the supercell case, which now
keeps the refined pass and the crystal's own lattice both in the configuration
that flipped and at HEAD; the axial-row arm still fires on the crystal it was
written for.

The post-refinement commit gate required each residual family to fall by any
amount. The positional family is what sees the detector geometry a commit hands
the next pass, and its sign was a coin toss where it did not move: on 9yzk it
read +7.1 % in one build and -0.14 % (0.1 of its paired noise) in the next, and
the second committed a geometry the positions had no opinion about (R-free
0.382 -> 0.412, radial misfit 0.30 -> 0.39). The gate now asks the positional
residual to fall by more than its paired standard error (both geometries read on
the same held-out reflections) and the excitation residual not to rise by more
than its own. On the 37-set validation this moves five fit decisions: 9yzk
(R-free 0.412 -> 0.394, radial misfit 0.39 -> 0.29, ISa 9.7 -> 7.8), two
canonical-pass probes whose output is unchanged, and one pre-pass commit on 9pbb
whose excitation moved +0.1 of its noise (ISa 14.77 -> 14.54, R-free 0.2353 ->
0.2357, inside the noise). The in-house myob, cytc, lyso and sparse sets are
md5-identical.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
2026-09-28 23:45:59 +02:00

332 lines
15 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 <cmath>
#include <filesystem>
#include "../common/DiffractionExperiment.h"
#include "../common/ScanResultGenerator.h"
#include "../writer/FileWriter.h"
#include "../reader/JFJochHDF5Reader.h"
#include "../rugnux/Rugnux.h"
#include "../rugnux/RugnuxCommandLine.h"
#include "../rugnux/SpotWidth.h"
namespace {
// Write a small VDS dataset of `n` flat images and return nothing (prefix_master.h5 +
// prefix_data_000001.h5 land in the test working directory).
void WriteTestDataset(const std::string &prefix, int n) {
RegisterHDF5Filter();
DiffractionExperiment x(DetJF(1));
x.FilePrefix(prefix).ImagesPerTrigger(n).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(n).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
x.BeamX_pxl(512).BeamY_pxl(256).DetectorDistance_mm(150).IncidentEnergy_keV(WVL_1A_IN_KEV)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
std::vector<int16_t> image(x.GetPixelsNum(), 5);
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
ScanResultGenerator generator(x);
for (int i = 0; i < n; i++) {
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
}
EndMessage end_message;
end_message.max_image_number = n;
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
}
TEST_CASE("Rugnux_AzInt", "[HDF5][Full]") {
WriteTestDataset("process_azint_in", 8);
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("process_azint_in_master.h5"));
auto dataset = reader.GetDataset();
REQUIRE(dataset);
ProcessConfig config;
config.mode = ProcessMode::AzimuthalIntegration;
config.nthreads = 2;
config.output_prefix = "process_azint_out";
Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config);
ProcessResult result;
REQUIRE_NOTHROW(result = process.Run());
CHECK_FALSE(result.cancelled);
CHECK(result.images_processed == 8);
REQUIRE(result.written_master_path.has_value());
{
// The _process.h5 links back to the source images and carries an azimuthal profile per image.
JFJochHDF5Reader out;
REQUIRE_NOTHROW(out.ReadFile("process_azint_out_process.h5"));
CHECK(out.GetNumberOfImages() == 8);
std::shared_ptr<JFJochReaderImage> img;
REQUIRE_NOTHROW(img = out.LoadImage(0));
REQUIRE(img);
CHECK_FALSE(img->ImageData().az_int_profile.empty());
}
reader.Close();
remove("process_azint_in_master.h5");
remove("process_azint_in_data_000001.h5");
remove("process_azint_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);
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("process_noout_in_master.h5"));
auto dataset = reader.GetDataset();
// Empty output prefix => process without writing any file.
ProcessConfig config;
config.mode = ProcessMode::AzimuthalIntegration;
config.nthreads = 3;
Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config);
auto result = process.Run();
CHECK_FALSE(result.cancelled);
CHECK(result.images_processed == 6);
CHECK_FALSE(result.written_master_path.has_value());
reader.Close();
remove("process_noout_in_master.h5");
remove("process_noout_in_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("Rugnux_Cancel", "[HDF5][Full]") {
WriteTestDataset("process_cancel_in", 8);
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("process_cancel_in_master.h5"));
auto dataset = reader.GetDataset();
ProcessConfig config;
config.mode = ProcessMode::AzimuthalIntegration;
config.nthreads = 2;
Rugnux process(reader, dataset->experiment, *dataset->pixel_mask, config);
process.Cancel(); // cancel before running: the worker loop stops immediately
auto result = process.Run();
CHECK(result.cancelled);
CHECK(result.images_processed == 0);
CHECK_FALSE(result.written_master_path.has_value());
reader.Close();
remove("process_cancel_in_master.h5");
remove("process_cancel_in_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("RugnuxCommandLine_Full", "[process]") {
DiffractionExperiment x(DetJF(1));
IndexingSettings idx;
idx.Algorithm(IndexingAlgorithmEnum::FFT);
idx.GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum::BeamCenter);
x.ImportIndexingSettings(idx);
x.SpaceGroupNumber(96);
ProcessConfig config;
config.mode = ProcessMode::FullAnalysis;
config.nthreads = 8;
config.output_prefix = "run1";
config.end_image = 500;
config.rotation_indexing = true;
config.two_pass_rotation = true;
config.rotation_indexing_image_count = 30;
config.spot_finding = DiffractionExperiment::DefaultDataProcessingSettings();
const std::string cmd = RugnuxCommandLine(config, x, "/data/test_master.h5");
CHECK(cmd.rfind("rugnux", 0) == 0);
CHECK(cmd.find("-N 8") != std::string::npos);
CHECK(cmd.find("-e 500") != std::string::npos);
CHECK(cmd.find("-o run1") != std::string::npos);
CHECK(cmd.find("-X fft") != std::string::npos);
CHECK(cmd.find("-S 96") != std::string::npos);
// -R takes an optional argument, so its value must be attached (-R30); a separate "-R 30" token
// would not re-parse (getopt would leave 30 as a positional and drop the count).
CHECK(cmd.find("-R30") != std::string::npos);
CHECK(cmd.find("-R 30") == std::string::npos);
CHECK(cmd.find("/data/test_master.h5") != std::string::npos);
}
TEST_CASE("RugnuxCommandLine_AzInt", "[process]") {
DiffractionExperiment x(DetJF(1));
AzimuthalIntegrationSettings a;
a.AzimuthalBinCount(4);
x.ImportAzimuthalIntegrationSettings(a);
ProcessConfig config;
config.mode = ProcessMode::AzimuthalIntegration;
config.nthreads = 2;
config.output_prefix = "az";
const std::string cmd = RugnuxCommandLine(config, x, "in.h5");
CHECK(cmd.rfind("rugnux", 0) == 0);
CHECK(cmd.find("--mode azint") != std::string::npos);
CHECK(cmd.find("--azim-phi-bins 4") != std::string::npos);
CHECK(cmd.find("--azim-min-q") != std::string::npos);
CHECK(cmd.find("in.h5") != std::string::npos);
}
namespace {
// A field of identical round Gaussian spots on three rings, so that the width estimator sees
// several resolution bands with the same true width and its 1/d fit has to come back flat.
void PaintGaussianSpots(ImagePreprocessorBuffer &image, int w, double sigma, double total_counts,
std::vector<DiffractionSpot> &spots) {
constexpr int BKG = 3;
for (size_t i = 0; i < image.size(); i++) image[i] = BKG;
const double amp = total_counts / (2.0 * M_PI * sigma * sigma);
for (int radius : {150, 350, 550})
for (int k = 0; k < 20; k++) {
const double phi = 2.0 * M_PI * k / 20.0 + 0.1 * radius;
const int cx = static_cast<int>(std::lround(600 + radius * std::cos(phi)));
const int cy = static_cast<int>(std::lround(600 + radius * std::sin(phi)));
for (int dy = -14; dy <= 14; dy++)
for (int dx = -14; dx <= 14; dx++)
image[static_cast<size_t>(cy + dy) * w + (cx + dx)] +=
static_cast<int32_t>(std::lround(
amp * std::exp(-(dx * dx + dy * dy) / (2.0 * sigma * sigma))));
spots.emplace_back(static_cast<uint32_t>(cx), static_cast<uint32_t>(cy),
static_cast<int64_t>(total_counts));
}
}
}
// The width the adaptive integration radius is set from. A round Gaussian of width sigma holds 80 %
// of its flux inside sqrt(2 ln 5) * sigma = 1.794 * sigma, and that is what the estimator has to
// return - over an aperture that owes nothing to the integrator's r1, which is the whole point of
// measuring it here rather than reading the integrator's own second moment.
TEST_CASE("SpotWidth_Gaussian", "[process]") {
constexpr int W = 1200, H = 1200;
DiffractionGeometry geometry;
geometry.BeamX_pxl(600).BeamY_pxl(600).DetectorDistance_mm(200).PixelSize_mm(0.075)
.Wavelength_A(1.0);
for (double sigma : {1.0, 2.2}) {
ImagePreprocessorBuffer image(static_cast<size_t>(W) * H);
std::vector<DiffractionSpot> spots;
PaintGaussianSpots(image, W, sigma, 20000.0, spots);
std::vector<spot_width::FluxCurve> curves;
MeasureSpotFluxCurves(image, W, H, geometry, spots, curves);
REQUIRE(curves.size() >= 45);
const auto r80 = spot_width::R80AtReference(curves);
REQUIRE(r80.has_value());
CHECK(*r80 == Catch::Approx(1.794 * sigma).margin(0.3));
}
// The rule the measurement drives: the shipped radius below the line, the capped one above it.
CHECK(spot_width::R1ForWidth(1.0f) == 4.0f);
CHECK(spot_width::R1ForWidth(1.794f) == 4.0f);
CHECK(spot_width::R1ForWidth(2.4f) == 5.0f);
CHECK(spot_width::R1ForWidth(3.947f) == 6.0f);
CHECK(spot_width::R1ForWidth(9.0f) == 6.0f);
}
// The rotation two-pass quality guard. The refined pass is sent back to the header geometry only for
// a decisive loss of signal - reflections at I/sigma >= 2 net of the noise tail at <= -2 - over no
// more of reciprocal space.
TEST_CASE("RefinedPassIsWorse", "[process]") {
const auto pass = [](const UnitCell &cell, char centering, int64_t strong, int64_t negative,
int64_t reflections) {
ProcessResult r;
r.has_merge_statistics = true;
r.consensus_cell = cell;
r.consensus_centering = centering;
r.search_merge_strong_reflections = strong;
r.search_merge_negative_reflections = negative;
r.search_merge_reflections = reflections;
return r;
};
const UnitCell cell{50.0f, 60.0f, 70.0f, 90.0f, 90.0f, 90.0f};
// Same lattice: a loss of more than 10 % of the signal is decisive, a smaller one is not.
const ProcessResult header = pass(cell, 'P', 50000, 1000, 200000);
CHECK(!RefinedPassIsWorse(header, pass(cell, 'P', 40000, 1000, 200000)).empty());
CHECK(RefinedPassIsWorse(header, pass(cell, 'P', 46000, 1000, 200000)).empty());
// ... unless the refined pass reaches into more of reciprocal space: a dilution, not a loss.
CHECK(RefinedPassIsWorse(header, pass(cell, 'P', 40000, 1000, 420000)).empty());
// The header pass on a 7x supercell: seven times the reflections, six in seven of them empty,
// and their noise tail alone reaches the strong count. Raw, 120k against 80k would read as a
// decisive loss for the refined pass on the crystal's own lattice; net of the noise it is none.
const UnitCell supercell{350.0f, 60.0f, 70.0f, 90.0f, 90.0f, 90.0f};
const ProcessResult on_supercell = pass(supercell, 'P', 120000, 40000, 1400000);
CHECK(RefinedPassIsWorse(on_supercell, pass(cell, 'P', 80000, 1500, 200000)).empty());
// The rescue the guard is there for is kept: a refined pass that drops to a sub-cell of the
// header's lattice loses the reflections the sub-cell cannot index. Half the volume, half the
// reflections - the same coverage - and half the signal.
const UnitCell subcell{25.0f, 60.0f, 70.0f, 90.0f, 90.0f, 90.0f};
CHECK(!RefinedPassIsWorse(header, pass(subcell, 'P', 25500, 500, 100000)).empty());
// A centred setting of the same lattice counts only the reflections its centring allows, so the
// counts compare as they are; the primitive volume puts the coverage on the same footing.
const UnitCell doubled{100.0f, 60.0f, 70.0f, 90.0f, 90.0f, 90.0f};
CHECK(RefinedPassIsWorse(header, pass(doubled, 'C', 49000, 1000, 200000)).empty());
CHECK(!RefinedPassIsWorse(header, pass(doubled, 'C', 40000, 1000, 200000)).empty());
// The axial arm: half of the low-order axial rows lost in the same setting, with no more signal.
ProcessResult lost_rows = pass(cell, 'P', 50000, 1000, 200000);
ProcessResult with_rows = header;
with_rows.search_merge_axial_reflections = 16;
lost_rows.search_merge_axial_reflections = 8;
CHECK(!RefinedPassIsWorse(with_rows, lost_rows).empty());
}
// The beam-centre arbitration between a first pass at the file's centre and one at the measured
// centre. Two arms on the same lattice are two geometries of one hypothesis: they are judged on the
// signal each measured, not on a CC1/2 that reads the same on both.
TEST_CASE("MeasuredCentreWins", "[process]") {
const auto arm = [](const UnitCell &cell, gemmi::CrystalSystem system, double cc_half, int64_t strong) {
ProcessResult r;
r.has_merge_statistics = true;
r.consensus_cell = cell;
r.consensus_centering = 'P';
r.rotation_lattice_type = LatticeMessage{'P', 0, system};
r.search_merge_cc_half = cc_half;
r.search_merge_strong_reflections = strong;
r.search_merge_reflections = 40000;
return r;
};
const UnitCell triclinic{40.0f, 41.0f, 100.0f, 86.0f, 84.0f, 72.0f};
// The same lattice in the setting with a and b exchanged, as a noisy a ~ b can come out.
const UnitCell swapped{41.02f, 39.98f, 100.1f, 84.0f, 86.0f, 72.0f};
const auto tri = gemmi::CrystalSystem::Triclinic;
const ProcessResult file = arm(triclinic, tri, 0.87, 9000);
REQUIRE(ArmsHoldSameLattice(file, arm(swapped, tri, 0.87, 9000)));
// Same lattice, same CC1/2: the arm that measured more signal wins, whichever centre it is...
CHECK(!MeasuredCentreWins(file, arm(swapped, tri, 0.87, 9500)).empty());
CHECK(MeasuredCentreWins(file, arm(swapped, tri, 0.87, 8500)).empty());
// ...a tie keeps the file's centre, and a higher CC1/2 alone does not move it.
CHECK(MeasuredCentreWins(file, arm(swapped, tri, 0.99, 9000)).empty());
// Different lattices: the metric question, decided on the search merges.
const UnitCell monoclinic{57.0f, 42.0f, 100.0f, 90.0f, 95.0f, 90.0f};
const auto mono = gemmi::CrystalSystem::Monoclinic;
REQUIRE(!ArmsHoldSameLattice(file, arm(monoclinic, mono, 0.87, 9000)));
CHECK(!MeasuredCentreWins(file, arm(monoclinic, mono, 0.95, 9000)).empty());
CHECK(MeasuredCentreWins(file, arm(monoclinic, mono, 0.88, 20000)).empty());
}
Reference in New Issue View Git Blame Copy Permalink