mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 88 Wiki Activity
Files
rc181
Jungfraujoch/tests/RugnuxTest.cpp
T
leonarski_fandClaude Opus 5.5 b822a775ac rugnux: solvent build optional, off by default (--refine-solvent / --no-refine-solvent) 
The refinement's solvent step (waters built and removed, ions relabelled, and the
solvent audit) now runs only with --refine-solvent; --no-refine-solvent is the
default. When off, no audit is made, no <prefix>_solvent.json is written, and the
report carries REFINE_SOLVENT=NOT_PERFORMED in place of the REFINE_SOLVENT_ keys.
RefinementSettings::solvent and RugnuxConfig::refine_solvent default to false, so
the broker's live rugnux follows; RefinementResult::solvent becomes optional.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqrcrEZy6XBUCMZvyT
2026-11-01 22:07:04 +01:00

972 lines
47 KiB
C++
Raw Permalink 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 "../reader/JFJochMultiReader.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_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);
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);
// A placed model is refined by default: only --no-refine is given back.
CHECK(RugnuxCommandLine(config, x, "in.h5", "", false).find("--no-refine") != std::string::npos);
config.refine = true;
CHECK(RugnuxCommandLine(config, x, "in.h5", "", false).find("--no-refine") == std::string::npos);
// The solvent build is off by default: only --refine-solvent is given back.
CHECK(RugnuxCommandLine(config, x, "in.h5", "", false).find("--refine-solvent") == std::string::npos);
config.refine_solvent = true;
CHECK(RugnuxCommandLine(config, x, "in.h5", "", false).find("--refine-solvent") != std::string::npos);
config.refine_solvent = false;
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());
// The weak-sweep arm speaks only where the probe measures the lattice in no shell.
CHECK(!DoubledCellWorthMergingWeak(breaks_down));
CHECK(!DoubledCellWorthMergingWeak(flat));
// Lattice unmeasured on the probe's frames: the class clears the null (z > 2 x 3) in two shells.
const auto weak = trend({{0.12, 2.4, 0.3}, {0.37, 2.8, 0.0}, {0.62, 6.8, 0.0}, {0.56, 7.0, 0.0}});
CHECK(!DoubledCellWorthMerging(weak));
CHECK(DoubledCellWorthMergingWeak(weak));
// ... and in one shell only: not enough.
const auto weak_one = trend({{0.12, 2.4, 0.3}, {0.37, 2.8, 0.0}, {0.62, 6.8, 0.0}, {0.56, 4.0, 0.0}});
CHECK(!DoubledCellWorthMergingWeak(weak_one));
}
TEST_CASE("MergedDoubledClassTrend", "[process]") {
// A doubled cell whose l-odd class is a fifth of the l-even intensity at low resolution and all of it
// at high: the trend reads the ratio per probe shell, and rises.
const UnitCell doubled{50, 60, 140, 90, 90, 90}, half{50, 60, 70, 90, 90, 90};
const auto parity = DoubledClassParity(doubled, half, 'P');
REQUIRE(parity.has_value());
CHECK((*parity == std::array<int, 3>{0, 0, 1}));
std::mt19937 rng(7);
std::exponential_distribution<double> wilson(1.0);
std::vector<MergedReflection> merged;
for (int h = 0; h <= 20; h++)
for (int k = 0; k <= 24; k++)
for (int l = 0; l <= 56; l++) {
const double s2 = h * h / 2500.0 + k * k / 3600.0 + l * l / 19600.0;
if (s2 <= 0.0 || s2 > 1.0 / (2.5 * 2.5))
continue;
MergedReflection r;
r.h = h; r.k = k; r.l = l;
r.d = static_cast<float>(1.0 / std::sqrt(s2));
const double t = std::clamp((1.0 / 8.0 - 1.0 / r.d) / (1.0 / 8.0 - 1.0 / 2.5), 0.0, 1.0);
const double share = (l % 2) ? 0.2 + 0.8 * t : 1.0;
r.I = static_cast<float>(100.0 * share * wilson(rng));
r.sigma = 5.0f;
merged.push_back(r);
}
ProcessResult::SupercellClass sc;
sc.trend = MergedDoubledClassTrend(merged, *parity);
REQUIRE(sc.trend.size() >= 4);
CHECK(sc.trend.front().excess < 0.4);
CHECK(sc.trend.back().excess > 0.7);
CHECK(!DoubledCellRises(sc).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);
// A small sample - a few hundred added reflections, below the probe arm's floor of 1000: refused at that
// floor, fitted without it, and then the share over all shells, two standard errors down, decides. Over
// all shells a flat 0.3 is too little to tell from noise on this many reflections; a class that is
// mostly flat (1.0 added) is told.
std::vector<MergedReflection> small_clean, small_contaminated;
for (size_t i = 0; i < clean.size(); i += 30) {
small_clean.push_back(clean[i]);
small_contaminated.push_back(clean[i]);
if (clean[i].l & 1)
small_contaminated.back().I += 1.0f;
}
CHECK(!FitDoubledClassTranslation(small_clean, doubled, half).tested);
const DoubledClassTranslation c = FitDoubledClassTranslation(small_clean, doubled, half, 'P', 0);
const DoubledClassTranslation e = FitDoubledClassTranslation(small_contaminated, doubled, half, 'P', 0);
REQUIRE(c.tested);
REQUIRE(e.tested);
CHECK(c.share_se > 0.0);
CHECK(c.share - 2.0 * c.share_se >= 0.5);
CHECK(e.share - 2.0 * e.share_se < 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));
}
TEST_CASE("ResolvableAxisLength", "[LongAxis]") {
// 500 mm, 1 A, 100 um pixels, 2 px spots: rows of a 2500 A axis are still a spot apart.
CHECK(ResolvableAxisLength_A(500.0f, 1.0f, 0.1f, 2.0f) == Catch::Approx(2500.0));
// A short distance and large pixels resolve only a protein-sized axis: the hypothesis is not asked.
CHECK(ResolvableAxisLength_A(180.0f, 1.0f, 0.172f, 3.0f) < 500.0f);
// Twice the distance, twice the axis; twice the spot width, half of it.
CHECK(ResolvableAxisLength_A(400.0f, 1.0f, 0.075f, 2.0f)
== Catch::Approx(2.0 * ResolvableAxisLength_A(200.0f, 1.0f, 0.075f, 2.0f)));
CHECK(ResolvableAxisLength_A(200.0f, 1.0f, 0.075f, 4.0f)
== Catch::Approx(0.5 * ResolvableAxisLength_A(200.0f, 1.0f, 0.075f, 2.0f)));
}
TEST_CASE("IsSuperLattice", "[LongAxis]") {
const CrystalLattice cell(60.0f, 70.0f, 200.0f, 90.0f, 95.0f, 90.0f);
const Coord a = cell.Vec0(), b = cell.Vec1(), c = cell.Vec2();
// An exact index-4 supercell in another basis: every node of it is a node of the cell.
const CrystalLattice super(a * 2.0f, b + a, c * 2.0f);
CHECK(IsSuperLattice(super, cell));
// The relation runs one way: the cell is a SUB-lattice of its supercell, not a superlattice.
CHECK_FALSE(IsSuperLattice(cell, super));
// The same lattice in another basis is not a superlattice of itself.
CHECK_FALSE(IsSuperLattice(CrystalLattice(a, b + a, c - b), cell));
// A cell 1.37x longer along c - the volume of a wrong long axis against the true one - is not
// related to it at all, in either direction.
const CrystalLattice longer(a, b, c * 1.37f);
CHECK_FALSE(IsSuperLattice(longer, cell));
CHECK_FALSE(IsSuperLattice(cell, longer));
// A 4.06x volume that is not exact (a triclinic cell near a 4x multiple) is not a superlattice.
CHECK_FALSE(IsSuperLattice(CrystalLattice(a * 2.0f, b + a * 0.5f, c * 2.03f), cell));
}
TEST_CASE("LongAxisChoice_PooledSpotsNotFrames", "[LongAxis]") {
// Shares of 10000 pooled spots; `reflections` is what the beats-chance test reads.
const auto ev = [](int64_t on, int64_t chance) {
return ValidationSpotEvidence{.spots = 10000, .on_lattice = on, .by_chance = chance,
.reflections = on, .reflections_by_chance = chance};
};
// The standing lattice takes 20% of the spots. A long-axis cell taking 60% is adopted.
const ValidationSpotEvidence standing = ev(2000, 60);
CHECK(LongAxisChoice(standing, {{ev(6000, 60), false}}) == std::optional<size_t>(0));
// Of two candidates the one with more spots beyond chance wins - not the first, not the longer.
CHECK(LongAxisChoice(standing, {{ev(3900, 50), false}, {ev(7800, 50), false}}) == std::optional<size_t>(1));
// One with a large share by chance is credited only the excess.
CHECK(LongAxisChoice(standing, {{ev(6000, 4500), false}, {ev(5000, 50), false}}) == std::optional<size_t>(1));
// Within the noise of the standing share, the standing lattice is kept.
CHECK_FALSE(LongAxisChoice(standing, {{ev(2050, 60), false}}).has_value());
// A superlattice of the standing cell is not adopted on the pooled count, which is biased toward it.
CHECK_FALSE(LongAxisChoice(standing, {{ev(6000, 60), true}}).has_value());
CHECK(LongAxisChoice(standing, {{ev(9000, 60), true}, {ev(5000, 60), false}}) == std::optional<size_t>(1));
// A standing pass with no lattice is zero evidence; a candidate must still beat chance.
CHECK(LongAxisChoice({}, {{ev(3000, 60), false}}) == std::optional<size_t>(0));
CHECK_FALSE(LongAxisChoice({}, {{ev(50, 60), false}}).has_value());
CHECK_FALSE(LongAxisChoice(standing, {}).has_value());
}
// The sub-lattice decision: the standing cell's spots on the nodes the sub-lattice drops, against what
// a wrong spindle angle puts on them.
TEST_CASE("SubLatticeSuffices", "[process]") {
// A sub-lattice that beats chance by far.
const ValidationSpotEvidence sub{.spots = 1400, .on_lattice = 340, .by_chance = 5, .reflections = 300,
.reflections_by_chance = 30};
// The dropped nodes hold what chance puts there, or a few more within its noise: halve.
CHECK(SubLatticeSuffices(3, 2.5, sub));
CHECK(SubLatticeSuffices(7, 2.5, sub)); // 4.5 over, 3.29 sqrt(2.5) = 5.2
// A class of its own (a superstructure, a pseudo-translation, or half of a correct protein cell's
// reflections): keep the standing cell.
CHECK_FALSE(SubLatticeSuffices(12, 2.5, sub));
CHECK_FALSE(SubLatticeSuffices(4000, 70.0, sub));
// A sub-lattice that does not beat chance is no hypothesis.
const ValidationSpotEvidence chance{.spots = 1400, .on_lattice = 10, .by_chance = 9, .reflections = 10,
.reflections_by_chance = 45};
CHECK_FALSE(SubLatticeSuffices(0, 2.5, chance));
}
// Two domains of a monoclinic cell twinned by a twofold about c*. With c cos(beta) = -a/4 the twin
// takes c to c + a/2, so the domains coincide on the reflections with h even and their coincidence
// lattice doubles c; any supercell of that one also indexes both. A synthetic cell.
TEST_CASE("TwinCompositeAt", "[process]") {
const float beta = std::acos(-1.0f / 6.0f) * 180.0f / static_cast<float>(M_PI);
const CrystalLattice domain(6.0f, 8.0f, 9.0f, 90.0f, beta, 90.0f);
const Coord a = domain.Vec0(), b = domain.Vec1(), c = domain.Vec2();
const Coord u = domain.Cstar() * (1.0f / domain.Cstar().Length());
const auto twin = [&](const Coord &q) { return u * (2.0f * (q * u)) - q; };
std::vector<Coord> cloud, single;
for (int h = -4; h <= 4; h++)
for (int k = -4; k <= 4; k++)
for (int l = -4; l <= 4; l++) {
if (h == 0 && k == 0 && l == 0)
continue;
const Coord q = domain.Astar() * h + domain.Bstar() * k + domain.Cstar() * l;
if (q.Length() > 1.0f)
continue;
cloud.push_back(q);
cloud.push_back(twin(q));
single.push_back(q);
}
const auto indices = [](const CrystalLattice &s, const std::vector<Coord> &qs) {
std::vector<std::array<int, 3>> hkl;
for (const auto &q : qs)
hkl.push_back({static_cast<int>(std::lround(s.Vec0() * q)), static_cast<int>(std::lround(s.Vec1() * q)),
static_cast<int>(std::lround(s.Vec2() * q))});
return hkl;
};
const auto same_cell = [](const CrystalLattice &x, const CrystalLattice &y) {
const auto ux = x.GetUnitCell(), uy = y.GetUnitCell();
return std::abs(ux.a - uy.a) < 0.01f && std::abs(ux.b - uy.b) < 0.01f && std::abs(ux.c - uy.c) < 0.01f
&& std::abs(std::abs(std::cos(ux.beta * M_PI / 180.0)) - std::abs(std::cos(uy.beta * M_PI / 180.0))) < 1e-3;
};
const CrystalLattice reduced = domain.NiggliReduce();
// The coincidence lattice, index 2.
const CrystalLattice doubled = CrystalLattice(a, b, c * 2.0f).NiggliReduce();
const auto at2 = TwinCompositeAt(doubled, indices(doubled, cloud), 2, 0.01);
REQUIRE(at2.has_value());
CHECK(at2->composite);
CHECK(at2->neither == 0);
CHECK(at2->p_both == Catch::Approx(0.25));
CHECK(same_cell(at2->domain, reduced));
// The twin law, a twofold up to the lattice's own rotations (about c*, or about a, which the
// monoclinic twofold about b takes it to): either reverses b, so 0k0 is its own partner's Friedel
// mate and is kept.
REQUIRE(at2->rotation.has_value());
const RotMatrix &r = *at2->rotation;
CHECK(r.Column(0).x + r.Column(1).y + r.Column(2).z == Catch::Approx(-1.0f).margin(1e-3));
CHECK(((r * b) + b).Length() < 1e-3f * b.Length());
// The twin fraction from the intensity of the spots on one domain alone: domain two at half the
// intensity of domain one is a third of the crystal (which domain comes first is the count's tie).
std::vector<float> intensity;
for (size_t i = 0; i < cloud.size(); i++)
intensity.push_back(i % 2 ? 0.5f : 1.0f);
const auto weighted = TwinCompositeAt(doubled, indices(doubled, cloud), 2, 0.01, intensity);
REQUIRE(weighted.has_value());
CHECK(std::min(weighted->fraction, 1.0 - weighted->fraction) == Catch::Approx(1.0 / 3.0).margin(1e-3));
// A skewed index-5 supercell of it: the domains are index-10 sub-lattices, among them forms with no
// coefficient prime to 10.
const CrystalLattice skewed = CrystalLattice(a * 2.0f + b, b * 2.0f - a, c * 2.0f).NiggliReduce();
const auto at10 = TwinCompositeAt(skewed, indices(skewed, cloud), 10, 0.01);
REQUIRE(at10.has_value());
CHECK(at10->composite);
CHECK(at10->neither == 0);
CHECK(same_cell(at10->domain, reduced));
// The same two domains with the class on neither holding far more than the chance share says it
// can (a share of 1% of the spots by chance): not two domains.
auto crowded = indices(doubled, cloud);
const size_t n_domains = crowded.size();
// Spots on nodes of every class, as strays fall: a quarter of the nodes are on neither.
std::mt19937 rng(7);
std::uniform_int_distribution<int> index(-6, 6);
for (int i = 0; i < 400; i++)
crowded.push_back({index(rng), index(rng), index(rng)});
const auto with_neither = TwinCompositeAt(doubled, crowded, 2, 0.01);
REQUIRE(with_neither.has_value());
CHECK(with_neither->neither_expected == Catch::Approx(0.01 * (n_domains + 400) * 0.25));
CHECK_FALSE(with_neither->composite);
// One domain alone on the doubled cell: one lattice, not a composite.
const auto alone = TwinCompositeAt(doubled, indices(doubled, single), 2, 0.01);
CHECK(!(alone.has_value() && alone->composite));
}
Reference in new issue View Git Blame Copy Permalink