mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 87 Wiki Activity
Files
0cf62452cd452b038e4e38fef7d0ac53258ab68f
Jungfraujoch/tests/RugnuxTest.cpp
T
leonarski_fandClaude Opus 5.5 d6166f6fa5 Panel layout stage 2: prediction and Bragg integration through the panel table 
Everything sparse that took a radial quantity from an image coordinate now takes it from the
reflection's (or spot's) position on the detector plane, through its panel (DetectorModules).
With an empty module table every changed expression reduces to the old one in the same
floating-point order, so flat detectors - every rotation set, EIGER, CCDs, assembled JUNGFRAU -
integrate bit for bit as before.

- Prediction (CPU, GPU, rotation CPU and GPU): bounds test against PlaneBounds instead of the
  image rectangle.
- IntegrateOnImage: a prediction in a gap or off every panel is dropped (PlaneToPanel), no longer
  passed through onto an unrelated stored pixel. The panel id stays local: the engine and the way
  back find it from the image position (the panel boxes tile the image), and a prediction that
  rounding to float put on another panel's edge is dropped so the two agree. Nothing is added to
  Reflection.
- Bragg engines (CPU, CUDA, Metal) take a BraggPanel per reflection (BraggStencil.h): its plane
  position, the panel's pixel axes on the plane and its box. The stencil's radial direction is
  taken on the plane and expressed in the panel's pixel axes; r0, the radial background curve,
  the Gaussian's elongation, the learned moments, tan(2theta) and Finalize's polarisation use the
  plane position; r_max / n_rad come from PlaneBounds. The neighbour mask and the owner map are
  clamped to the reflection's own panel; inside the walked box a pixel off the reflection's panel
  is read as a gap (unreadable, statically masked - not an overload), which is what the assembled
  image has there, so a disk cut by a panel edge is partial rather than silently "full".
- Risk R4: the ProfileEmpirical grid is learned and applied in the plane's axes (the panel axes
  taken to the nearest quarter turn, BraggProfileCell), so turned panels pool one shape. Without
  it the stacked layout read up to 0.3 % off the assembled one in the harness.
- SpotWidth / SpotFootprint: aperture bounds against the spot's panel box, beam distance,
  resolution and the shape's 2theta Jacobian on the plane, moments read in the panel's axes.
- GetDetectorMaxQ_recipA, the far-corner radii (first-pass tilt report, metric-symmetry arm,
  footprint table range) from PlaneBounds; pre-scan spot resolution on the plane;
  LogHeaderBeamCenter: geometric centre of the panels' plane box, and the masked-pixel check on
  the panel the beam lands on (none when it is in a gap).

Harness: PanelLayout_BraggIntegration is now a real [PanelLayout] check, over six settings
(default, bandwidth-elongated ring, radial background, overlap exclude, box sum, empirical
profile, all with polarisation), stacked vs assembled to 1e-5; PanelLayout_BraggIntegrationGPU
runs the same on the CUDA engine. Its spots avoid exact half pixels: the grid centres on the
rounded position, and a half rounds to the other neighbour on a turned panel.

The Metal engine (.metal/.mm) mirrors the CPU change exactly but is UNCOMPILED here (no macOS):
it needs a build and BraggIntegrationEngineMetal_BitIdenticalToCPU on a Mac.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
2026-10-21 10:40:20 +02:00

693 lines
32 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 <random>
#include <catch2/catch_all.hpp>
#include <cmath>
#include <filesystem>
#include <mutex>
#include <set>
#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/RugnuxOutputFiles.h"
#include "../rugnux/RugnuxDefaults.h"
#include "../image_analysis/bragg_integration/SpotWidth.h"
#include "../image_analysis/geom_refinement/Calibrants.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, bool rotation = false) {
RegisterHDF5Filter();
DiffractionExperiment x(DetJF(1));
if (rotation)
x.Goniometer(GoniometerAxis("omega", 0.0f, 1.0f, Coord(1, 0, 0), std::nullopt));
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;
// The per-image analysis goes to the caller's sink, which writes it.
ProcessFileSink sink("process_azint_out");
config.image_sink = &sink;
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);
{
// 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();
// No image sink => nothing is streamed anywhere.
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);
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);
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);
}
// A rotation run that finds no lattice still reports every image - its spots, resolution and background,
// and that it did not index - as the canonical pass, and hands the same back with the failure.
TEST_CASE("Rugnux_NoLatticeReportsEveryImage", "[HDF5][Full]") {
WriteTestDataset("process_nolattice_in", 8, true);
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("process_nolattice_in_master.h5"));
auto dataset = reader.GetDataset();
DiffractionExperiment x = dataset->experiment;
auto indexing = x.GetIndexingSettings();
indexing.RotationIndexing(true);
indexing.Algorithm(IndexingAlgorithmEnum::Auto);
x.ImportIndexingSettings(indexing);
REQUIRE(x.IsRotationIndexing());
struct Observer : RugnuxObserver {
std::atomic<bool> canonical{false};
std::mutex m;
std::set<int64_t> images;
bool complete = true;
void OnImagePass(bool c) override { canonical = c; }
void OnImageProcessed(const DataMessage &msg) override {
if (!canonical)
return;
std::lock_guard lock(m);
images.insert(msg.number);
complete &= msg.spot_count.has_value() && msg.bkg_estimate.has_value()
&& msg.indexing_result == false;
}
} observer;
ProcessConfig config;
config.nthreads = 2;
config.rotation_indexing = true;
config.rotation_postrefine_geometry = true;
config.spot_finding = RugnuxDefaultSpotFinding(true);
Rugnux process(reader, x, *dataset->pixel_mask, config);
std::shared_ptr<ProcessResult> per_image;
try {
process.Run(&observer);
} catch (const NoLatticeFound &e) {
per_image = e.per_image;
}
REQUIRE(per_image);
CHECK(per_image->images_processed == 8);
CHECK(per_image->end_message.max_image_number == 8);
REQUIRE(per_image->outputs);
CHECK(per_image->outputs->process_start.has_value());
CHECK(per_image->end_message.spot_count.size() == 8);
CHECK(per_image->end_message.v_bkg_estimate.size() == 8);
CHECK(observer.images.size() == 8);
CHECK(observer.complete);
reader.Close();
remove("process_nolattice_in_master.h5");
remove("process_nolattice_in_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// A calibration re-integrates its images binned about the geometry it fitted, and fits again. That
// second profile must be the sum over every image: here each of six frames carries a different sixth
// of the LaB6 rings, so a profile missing a frame is missing rings. The header's beam centre is 5 px
// out, far enough that the first pass's smeared profile loses to the re-binned one - so the geometry
// that comes out is the one fitted on the re-binned sum, and it must not depend on the thread count.
TEST_CASE("Rugnux_CalibrationRebinsEveryImage", "[HDF5][DetGeomCalib]") {
RegisterHDF5Filter();
constexpr int n = 6;
constexpr float true_beam_x = 1090.0f, true_beam_y = 1100.0f;
const std::vector<float> rings_q = CalibrantRings("lab6");
DiffractionExperiment x(DetJF4M());
x.FilePrefix("process_calib_in").ImagesPerTrigger(n).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(n).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
x.BeamX_pxl(true_beam_x).BeamY_pxl(true_beam_y).DetectorDistance_mm(100).IncidentEnergy_keV(WVL_1A_IN_KEV)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
const DiffractionGeometry geom_true = x.GetDiffractionGeometry();
x.BeamX_pxl(true_beam_x + 5.0f); // what the header says
// Radii of the rings the beam reaches, in pixels about the true centre.
std::vector<float> ring_radius;
for (const float q : rings_q) {
const float d = static_cast<float>(2.0 * PI) / q;
if (d > geom_true.GetWavelength_A() / 2.0f)
ring_radius.push_back(geom_true.ResToPxl(d));
}
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
ScanResultGenerator generator(x);
std::vector<int16_t> image(x.GetPixelsNum());
for (int i = 0; i < n; i++) {
for (int64_t p = 0; p < x.GetPixelsNum(); p++) {
const float px = static_cast<float>(p % x.GetXPixelsNum());
const float py = static_cast<float>(p / x.GetXPixelsNum());
const float r = std::hypot(px - true_beam_x, py - true_beam_y);
float v = 5.0f;
for (size_t k = i; k < ring_radius.size(); k += n) {
const float dr = r - ring_radius[k];
v += 400.0f * std::exp(-0.5f * dr * dr / (1.5f * 1.5f));
}
image[p] = static_cast<int16_t>(std::lround(v));
}
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();
}
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("process_calib_in_master.h5"));
auto dataset = reader.GetDataset();
REQUIRE(dataset);
DiffractionExperiment experiment(dataset->experiment);
auto azint = experiment.GetAzimuthalIntegrationSettings();
azint.AzimuthalBinCount(32);
experiment.ImportAzimuthalIntegrationSettings(azint);
std::vector<CalibrationResult> results;
for (const int nthreads : {1, 4}) {
ProcessConfig config;
config.mode = ProcessMode::Calibration;
config.calibrant_ring_q = rings_q;
config.nthreads = nthreads;
config.beam_center_check = false;
config.spot_finding = DiffractionExperiment::DefaultDataProcessingSettings();
config.spot_finding.indexing = false;
Rugnux process(reader, experiment, *dataset->pixel_mask, config);
ProcessResult result;
REQUIRE_NOTHROW(result = process.Run());
REQUIRE(result.calibration.has_value());
results.push_back(*result.calibration);
}
for (const auto &cal : results) {
CHECK(cal.converged);
CHECK(cal.geometry.GetBeamX_pxl() == Catch::Approx(true_beam_x).margin(0.5));
CHECK(cal.geometry.GetBeamY_pxl() == Catch::Approx(true_beam_y).margin(0.5));
}
CHECK(results[0].ring_points == results[1].ring_points);
CHECK(results[0].rms_radial_pxl == results[1].rms_radial_pxl);
CHECK(results[0].geometry.GetBeamX_pxl() == results[1].geometry.GetBeamX_pxl());
CHECK(results[0].geometry.GetBeamY_pxl() == results[1].geometry.GetBeamY_pxl());
CHECK(results[0].geometry.GetDetectorDistance_mm() == results[1].geometry.GetDetectorDistance_mm());
reader.Close();
remove("process_calib_in_master.h5");
remove("process_calib_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.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", "run1", false);
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 cufft") != 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);
CHECK(cmd.find("--model") == std::string::npos);
// --mr is --model with the search forced, and is given back as itself.
auto model = std::make_shared<RugnuxModel>();
model->name = "model.cif";
config.model = model;
CHECK(RugnuxCommandLine(config, x, "in.h5", "", false).find("--model model.cif") != std::string::npos);
config.force_mr = true;
const std::string mr_cmd = RugnuxCommandLine(config, x, "in.h5", "", false);
CHECK(mr_cmd.find("--mr model.cif") != std::string::npos);
CHECK(mr_cmd.find("--model") == std::string::npos);
}
TEST_CASE("RugnuxDefaults_BareRun", "[process]") {
// What both front ends (the rugnux CLI and the viewer's "Analyze dataset") start a run from.
const SpotFindingSettings rot = RugnuxDefaultSpotFinding(true);
const SpotFindingSettings stills = RugnuxDefaultSpotFinding(false);
CHECK(rot.adaptive_threshold);
CHECK(stills.adaptive_threshold);
CHECK(rot.min_pix_per_spot == 2);
CHECK(!stills.min_pix_per_spot.has_value());
CHECK(!rot.high_resolution_limit.has_value());
CHECK(rot.indexing);
DiffractionExperiment x(DetJF(1));
ApplyRugnuxExperimentDefaults(x);
CHECK(x.GetPolarizationFactor() == 0.99f);
CHECK(!x.GetBraggIntegrationSettings().GetBackgroundRadialCorrection().has_value()); // auto
DiffractionExperiment rot_x = x, stills_x = x;
ApplyRugnuxIntegrationRadii(rot_x, true);
ApplyRugnuxIntegrationRadii(stills_x, false);
CHECK(rot_x.GetBraggIntegrationSettings().GetR1() == BraggIntegrationSettings().GetR1());
CHECK(stills_x.GetBraggIntegrationSettings().GetR1() == 6.0f);
CHECK(stills_x.GetBraggIntegrationSettings().GetR3() == 14.0f);
CHECK(RugnuxDetectIceRings(std::nullopt, true));
CHECK(!RugnuxDetectIceRings(std::nullopt, false));
CHECK(!RugnuxDetectIceRings(false, true));
// A stationary axis is stills, the same test everywhere the workflow is decided.
DiffractionExperiment still_axis = x, turning_axis = x;
still_axis.Goniometer(GoniometerAxis("omega", 0.0f, 0.0f, Coord(1, 0, 0), std::nullopt));
turning_axis.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(1, 0, 0), std::nullopt));
CHECK(!RugnuxRotationData(x));
CHECK(!RugnuxRotationData(still_axis));
CHECK(RugnuxRotationData(turning_axis));
}
TEST_CASE("RugnuxCommandLine_AutomaticChoices", "[process]") {
// A run left to rugnux's automatic choices must give back a command line that leaves them automatic:
// -N, --max-spots and --integration-radius each switch one of them off.
DiffractionExperiment x(DetJF(1));
ProcessConfig config;
config.mode = ProcessMode::FullAnalysis;
config.nthreads = 8;
config.nthreads_auto = true;
config.measure_spot_budget = true;
config.adaptive_integration_radius = true;
config.rotation_indexing = true;
const std::string cmd = RugnuxCommandLine(config, x, "in.h5", "", false);
CHECK(cmd.find("-N ") == std::string::npos);
CHECK(cmd.find("--max-spots") == std::string::npos);
CHECK(cmd.find("--integration-radius") == std::string::npos);
config.nthreads_auto = false;
config.measure_spot_budget = false;
config.adaptive_integration_radius = false;
const std::string pinned = RugnuxCommandLine(config, x, "in.h5", "", false);
CHECK(pinned.find("-N 8") != std::string::npos);
CHECK(pinned.find("--max-spots") != std::string::npos);
CHECK(pinned.find("--integration-radius") != 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;
const std::string cmd = RugnuxCommandLine(config, x, "in.h5", "az", false);
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, DetectorModules(), 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("DoubledCellDecision", "[process]") {
// Shells of the supercell probe's trend: d range, excess over the null, z, the strongest null z and
// the lattice's <I/sigma>.
const auto trend = [](std::vector<std::array<double, 3>> excess_z_lattice) {
ProcessResult::SupercellClass sc;
float d = 20.0f;
for (const auto &[excess, z, lattice] : excess_z_lattice) {
ProcessResult::SupercellClass::Shell sh;
sh.d_low_A = d;
sh.d_high_A = d *= 0.7f;
sh.excess = excess;
sh.z = z;
sh.sibling_z = 3.0;
sh.lattice_i_over_sigma = lattice;
sh.resolved = 0.9;
sc.trend.push_back(sh);
}
return sc;
};
// A pseudo-translation that breaks down: weak at low resolution, a third of the lattice where the
// lattice is still measured.
const auto breaks_down = trend({{0.02, 17, 18}, {0.10, 38, 7}, {0.26, 36, 6}, {0.35, 45, 2.4}, {0.46, 33, 0.8}});
CHECK(DoubledCellWorthMerging(breaks_down));
CHECK(!DoubledCellRises(breaks_down).empty());
// A class that stays a tenth of the lattice at every resolution: merged as well, but it does not rise.
const auto flat = trend({{0.09, 10, 70}, {0.11, 23, 28}, {0.10, 24, 33}, {0.11, 24, 18}, {0.10, 23, 8}});
CHECK(DoubledCellWorthMerging(flat));
CHECK(DoubledCellRises(flat).empty());
// A class no stronger than the null: nothing.
const auto null_class = trend({{0.01, 4, 20}, {0.02, 5, 10}, {0.03, 2, 5}, {0.05, 1, 2}});
CHECK(!DoubledCellWorthMerging(null_class));
// A class that rises significantly but stays a few percent of the lattice: the half cell holds.
const auto small = trend({{0.009, 40, 500}, {0.014, 60, 300}, {0.016, 90, 170}, {0.03, 60, 50}, {0.06, 40, 2}});
CHECK(DoubledCellRises(small).empty());
}
TEST_CASE("FitDoubledClassTranslation", "[process]") {
// A doubled orthorhombic cell whose two halves are one copy translated by c/2 + (0.03, 0, 0): the
// reflections with l odd carry |F|^2 sin^2(pi h.eps), the others cos^2. With nothing else on the
// added class the fit finds the translation and leaves nothing over; with a flat 0.3 added to
// the class, most of the class is left unexplained.
const UnitCell doubled{50, 60, 140, 90, 90, 90}, half{50, 60, 70, 90, 90, 90};
std::mt19937 rng(1);
std::exponential_distribution<double> wilson(1.0);
std::vector<MergedReflection> clean, contaminated;
for (int h = -25; h <= 25; h++)
for (int k = 0; k <= 30; k++)
for (int l = 0; l <= 70; l++) {
if (h == 0 && k == 0 && l == 0)
continue;
MergedReflection r;
r.h = h; r.k = k; r.l = l;
r.d = 1.0f / std::sqrt(static_cast<float>(h * h / 2500.0 + k * k / 3600.0 + l * l / 19600.0));
if (r.d < 2.0f)
continue;
const double base = wilson(rng), x = std::pow(std::sin(PI * 0.03 * h), 2);
r.I = static_cast<float>(base * ((l & 1) ? x : 1.0 - x));
clean.push_back(r);
if (l & 1)
r.I += 0.3f;
contaminated.push_back(r);
}
const DoubledClassTranslation a = FitDoubledClassTranslation(clean, doubled, half);
REQUIRE(a.tested);
CHECK(a.parity == std::array<int, 3>{0, 0, 1});
const DoubledClassTranslation b = FitDoubledClassTranslation(contaminated, doubled, half);
REQUIRE(b.tested);
const auto weakest = [](const DoubledClassTranslation &t) {
double w = INFINITY;
for (const auto &sh : t.shells)
w = std::min(w, sh.modulated_share);
return w;
};
CHECK(weakest(a) > 0.8);
CHECK(weakest(b) < 0.5);
}
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());
// A refined pass that left a twin overlap of index 3 out of integration merged two thirds of the
// nodes: two thirds of the signal is no loss, and its axial rows are not the header pass's rows.
ProcessResult twin = pass(cell, 'P', 34000, 700, 132000);
twin.twin_overlap_index = 3;
twin.search_merge_axial_reflections = 10;
CHECK(RefinedPassIsWorse(with_rows, twin).empty());
twin.search_merge_strong_reflections = 28000;
CHECK(!RefinedPassIsWorse(with_rows, twin).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());
// ...but not by less than twice the Poisson noise of the two counts (sqrt(18100) = 135)...
CHECK(MeasuredCentreWins(file, arm(swapped, tri, 0.87, 9100)).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());
}
TEST_CASE("MeasuredCentreClearlyLoses", "[process]") {
// Same lattice: the measured centre loses only where it is behind on the signal by more than the
// Poisson noise of the two counts AND its inner shell does not agree better by the margin.
const SearchMergeSignal file{.strong = 9000, .reflections = 40000, .cc_half = 0.9,
.inner_cc_half = 0.98, .inner_reflections = 2000};
SearchMergeSignal m = file;
m.strong = 8500; // 500 behind, noise sqrt(17500) = 132
CHECK(MeasuredCentreClearlyLoses(file, m, true, 2.0));
m.strong = 8800; // 200 behind: within two sigmas
CHECK_FALSE(MeasuredCentreClearlyLoses(file, m, true, 2.0));
// A file centre whose few inner reflections disagree: the measured centre's inner shell might
// still carry it, so behind on the signal is not yet a loss.
SearchMergeSignal inconsistent = file;
inconsistent.inner_cc_half = 0.6;
inconsistent.inner_reflections = 20;
m.strong = 8500;
m.inner_cc_half = 0.64;
m.inner_reflections = 20;
CHECK_FALSE(MeasuredCentreClearlyLoses(inconsistent, m, true, 2.0));
// Different lattices: the count alone - the subset's CC1/2 is not read.
SearchMergeSignal d = file;
d.strong = 3000;
d.cc_half = 0.99;
CHECK(MeasuredCentreClearlyLoses(file, d, false, 2.0));
d.strong = 8900; // 100 behind: within two sigmas
d.cc_half = 0.50;
CHECK_FALSE(MeasuredCentreClearlyLoses(file, d, false, 2.0));
// A merge with no strong reflection in it measured none of the crystal. (An arm without a merge
// is never asked: RunAllPasses judges it on all the frames.)
CHECK(MeasuredCentreClearlyLoses(file, SearchMergeSignal{}, false, 2.0));
}
Reference in new issue View Git Blame Copy Permalink