mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 89 Wiki Activity
Files
rc182
Jungfraujoch/tests/RugnuxTest.cpp
T
leonarski_fandClaude Opus 5.5 7907d2980e rugnux: decide a cell against its integer multiple on integrated cosets 
A cell and its axis multiple were decided by spot counts in several places,
which cannot work: the larger lattice contains the smaller, so it fits every
spot at least as well. The decision now sits in one place, on the integrated
intensity of the reflections the larger cell adds.

- First pass: a pair of scheme answers whose primitive volumes are 2-4x apart
  and that tie on validation frames keeps the SMALLER cell. The spot-count
  arbiter (larger vs smaller and vs own sub-lattice, with its late re-check)
  is gone; the index-2 larger cell is recorded as an offered hypothesis.
- First-pass ladder: the spot-count harmonic refusal of a winning rung is gone.
- Pass-1 supercell guard and the pass-2 supercell-collapse guard lose their
  "the arbiter settled the larger cell" exemption (larger_settled,
  supercell_settled_), which no longer exists.
- Supercell probe: every half-integer class is read against its own sibling
  null, all classes that stand clear are listed, and group closure of the two
  strongest (their sum must stand clear too) is reported.
- Doubled-cell hypothesis: adopts the doubled cell also where the class is part
  of the lattice (a quarter of the lattice's intensity or more in every
  measured shell, coarsest included), and, for a doubled cell an indexing
  scheme offered, where the class grows to a quarter of the lattice's
  intensity, without the one-pseudo-translation test. A superstructure only the
  probe has seen is decided exactly as before.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu
2026-11-03 12:08:44 +01:00

981 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());
// A class that is part of the lattice: as strong as the lattice in every shell, flat. It does not rise,
// and the half cell misses it at every resolution.
const auto lattice = trend({{0.95, 40, 30}, {1.05, 60, 20}, {0.90, 50, 9}, {1.00, 30, 3}});
CHECK(DoubledCellWorthMerging(lattice));
CHECK(DoubledCellRises(lattice).empty());
CHECK(!DoubledCellIsLattice(lattice).empty());
// A superstructure that breaks down is not that (weak at low resolution), nor is a flat tenth.
CHECK(DoubledCellIsLattice(breaks_down).empty());
CHECK(DoubledCellIsLattice(flat).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