mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 74 Wiki Activity
Files
a4be5e6cd0e4bdfec1fd83e5f885d7a9d9f9b4b1
Jungfraujoch/tests/ResultReportTest.cpp
T
leonarski_fandClaude Opus 5 a4be5e6cd0 diagnostics: a translational pseudo-symmetry is detected and named 
A pseudo-translation makes one class of reflections systematically strong and the
complementary class weak. It is the classic predictor of a failed molecular
replacement, it raises the second moment where twinning lowers it so each masks the
other's test, and rugnux did not look for it at all.

Look for it: a native Patterson from the merged intensities, the largest off-origin
peak taken as a fraction of the origin, and the vector refined against the data.
Both halves are scored against a null computed for that crystal rather than a bound
taken from elsewhere - the peak against intensities permuted within resolution
shells, the modulation depth against the same search started from random vectors -
because the noise floor of the first runs from 1.35 to 17.7 per cent across the
corpus, so no fixed threshold could mean the same thing twice. Requiring both is
what keeps the false-positive rate down: either alone calls one crystal in ten.

A translation the data are exactly invariant under is a lattice vector, not a
pseudo-symmetry. Recognising that removes eleven false calls from ninety-five merges
folded in P1 - each reading 83 to 102 per cent of origin against a genuine maximum of
62 - and adds one true positive whose pseudo-translation lies underneath its own
centring and is unreachable without it.

Report-only: it gates nothing and changes no reflection, no scale and no group. It
fires on eight of a hundred and thirty-seven crystals, and on those the verdict moves
from OK to WARNINGS, which is the point of it. Median cost 55 ms, worst 1.9 s on the
largest merge here.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
2026-09-05 20:08:40 +02:00

610 lines
31 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include "../common/DiffractionExperiment.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../reader/JFJochHDF5Reader.h"
#include "../rugnux/ResultReport.h"
#include "../common/CUDAWrapper.h"
#include "../writer/FileWriter.h"
// Two synthetic ranges, one of each shape the report has to handle.
namespace {
SweepQuality TestSweepQuality() {
SweepQuality out;
out.measured = true;
out.sweep_deg = 60.0f;
out.ranges.push_back(SweepQualityRange{.first_image = 100, .last_image = 149,
.reason = SweepQualityReason::CrystalOutOfBeam,
.severity = 0.83f, .rotation_deg = 5.0f,
.mean_relative_scale = 0.12f, .mean_relative_cc = 0.30f,
.indexed_fraction = 0.015f});
out.ranges.push_back(SweepQualityRange{.first_image = 400, .last_image = 499,
.reason = SweepQualityReason::RadiationDamage,
.severity = 0.41f, .rotation_deg = 10.0f,
.mean_relative_scale = 0.55f, .mean_relative_cc = 0.80f,
.indexed_fraction = 0.62f});
return out;
}
std::vector<std::string> ReasonVocabulary() {
std::vector<std::string> out;
for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
out.emplace_back(SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
return out;
}
}
TEST_CASE("SweepQuality_ReasonVocabulary", "[Diagnostics]") {
// The codes are an interface - they are written verbatim into <prefix>_report.txt and into HDF5,
// where the per-image codes index this list from 1. A rename or a reorder silently breaks every
// consumer, so pin both the spelling and the order here.
const auto codes = ReasonVocabulary();
REQUIRE(codes.size() == 5);
CHECK(codes[0] == "no_diffraction");
CHECK(codes[1] == "crystal_out_of_beam");
CHECK(codes[2] == "weak_diffraction");
CHECK(codes[3] == "loss_of_centring");
CHECK(codes[4] == "radiation_damage");
}
TEST_CASE("ResultReport_Render", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(600);
ProcessResult result;
result.images_processed = 600;
result.indexing_rate = 0.87f;
result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
result.space_group = *gemmi::find_spacegroup_by_number(96);
result.used_beam_x_pxl = 766.62f;
result.used_beam_y_pxl = 846.87f;
result.used_distance_mm = 243.53f;
result.has_merge_statistics = true;
result.merge_statistics.sweep_quality = TestSweepQuality();
const auto text = RenderResultReport("prefix", "in.h5", x, result);
// The stable keys a consumer greps for. The version is pinned on purpose: a key added to the
// report is a contract change, and this line is where it has to be acknowledged.
CHECK(text.find("\nREPORT_VERSION= 10\n") != std::string::npos);
// The flight-path medium is an ASSUMPTION - no file states it - so the report has to say which
// one was made and what it was worth, or the user cannot know a choice was taken on their behalf.
CHECK(text.find("\nFLIGHT_PATH= AIR\n") != std::string::npos);
CHECK(text.find("\nFLIGHT_PATH_WILSON_B= ") != std::string::npos);
CHECK(text.find("\nOUTPUT_PREFIX= prefix\n") != std::string::npos);
CHECK(text.find("\nINDEXING_RATE= 0.8700\n") != std::string::npos);
CHECK(text.find("\nSPACE_GROUP_NUMBER= 96\n") != std::string::npos);
// What the data could NOT decide, beside what they did. The group here was given rather
// than searched for, and it is one of the 22 that come in enantiomorphic pairs.
CHECK(text.find("\nSPACE_GROUP_ALTERNATIVES= NONE\n") != std::string::npos);
CHECK(text.find("\nSPACE_GROUP_ENANTIOMORPH= GIVEN\n") != std::string::npos);
// Written only where there WAS a refusal: "NONE" on every other run answers a question nobody
// asked. The vocabulary and the status live in the developer report now, and are pinned there.
CHECK(text.find("SPACE_GROUP_REFUSED_POINT_GROUP=") == std::string::npos);
CHECK(text.find("\nSWEEP_QUALITY_COUNT= 2\n") != std::string::npos);
CHECK(text.find("SWEEP_QUALITY_STATUS=") == std::string::npos);
CHECK(text.find("SWEEP_QUALITY_REASONS=") == std::string::npos);
RunProvenance dev;
dev.developer = true;
const auto dev_text = RenderResultReport("prefix", "in.h5", x, result, dev);
CHECK(dev_text.find("\nSWEEP_QUALITY_STATUS= COMPUTED\n") != std::string::npos);
CHECK(dev_text.find("\nSWEEP_QUALITY_REASONS= no_diffraction crystal_out_of_beam weak_diffraction "
"loss_of_centring radiation_damage\n") != std::string::npos);
// One table row per range, with the reason code verbatim.
CHECK(text.find(" 100 149 50 5.0 crystal_out_of_beam ")
!= std::string::npos);
CHECK(text.find(" 400 499 100 10.0 radiation_damage ")
!= std::string::npos);
// ... and one plain-English WARNING line per range, greppable by the marker alone. They sit in
// the SUMMARY at the top of the file now, which is what the verdict is composed from.
CHECK(text.find("\nWARNING_COUNT= 2\n") != std::string::npos);
CHECK(text.find("\nVERDICT= WARNINGS\n") != std::string::npos);
CHECK(text.find("\nPATHOLOGY_FLAGS= SWEEP_GAPS\n") != std::string::npos);
CHECK(text.find("VERDICT=") < text.find("SWEEP_QUALITY_COUNT="));
CHECK(text.find("\nWARNING: Frames 100-149 out of beam (5.0 deg,") != std::string::npos);
CHECK(text.find("\nWARNING: Frames 400-499 radiation damage (10.0 deg,") != std::string::npos);
}
TEST_CASE("ResultReport_RenderEmpty", "[Diagnostics]") {
// A clean run and a run that never looked must be distinguishable: both have a count of 0, and
// only the STATUS key separates them. This is the property a consumer relies on.
DiffractionExperiment x(DetJF(1));
ProcessResult clean;
clean.has_merge_statistics = true;
clean.merge_statistics.sweep_quality.measured = true;
RunProvenance dev;
dev.developer = true;
const auto clean_text = RenderResultReport("p", "in.h5", x, clean);
CHECK(clean_text.find("\nSWEEP_QUALITY_COUNT= 0\n") != std::string::npos);
CHECK(RenderResultReport("p", "in.h5", x, clean, dev)
.find("\nSWEEP_QUALITY_STATUS= COMPUTED\n") != std::string::npos);
// An empty table used to print its header and two rules around nothing on every clean run.
CHECK(clean_text.find("FIRST_IMAGE LAST_IMAGE") == std::string::npos);
CHECK(clean_text.find("no degraded ranges") != std::string::npos);
ProcessResult not_merged;
const auto not_merged_text = RenderResultReport("p", "in.h5", x, not_merged);
CHECK(RenderResultReport("p", "in.h5", x, not_merged, dev)
.find("\nSWEEP_QUALITY_STATUS= NOT_COMPUTED\n") != std::string::npos);
CHECK(not_merged_text.find("\nMERGE= NOT_PERFORMED\n") != std::string::npos);
CHECK(not_merged_text.find("FIRST_IMAGE LAST_IMAGE") == std::string::npos);
// No lattice was found, and that is the FIRST thing the report says.
CHECK(not_merged_text.find("\nVERDICT= FAILED\n") != std::string::npos);
}
TEST_CASE("ResultReport_RadiationDamage", "[Diagnostics]") {
// RADIATION_DAMAGE_RELATIVE_B is a number only when there is one. A curve that was measured but
// that no straight line describes, and a monitor that could not run at all, are different answers,
// and a consumer has to be able to tell them apart - and both from a measured zero.
DiffractionExperiment x(DetJF(1));
ProcessResult result;
result.has_merge_statistics = true;
result.radiation_damage_text = "per-batch relative-B";
result.merge_statistics.radiation_damage_batch_deg = 10.0;
result.merge_statistics.radiation_damage_b_batch = {0.0f, 4.0f, NAN};
result.merge_statistics.radiation_damage_delta_b = 8.25;
CHECK(RenderResultReport("p", "in.h5", x, result).find("\nRADIATION_DAMAGE_RELATIVE_B= 8.25\n")
!= std::string::npos);
result.merge_statistics.radiation_damage_delta_b = NAN;
CHECK(RenderResultReport("p", "in.h5", x, result).find("\nRADIATION_DAMAGE_RELATIVE_B= NOT_A_TREND\n")
!= std::string::npos);
result.merge_statistics.radiation_damage_b_batch.clear();
CHECK(RenderResultReport("p", "in.h5", x, result).find("\nRADIATION_DAMAGE_RELATIVE_B= NOT_MEASURED\n")
!= std::string::npos);
}
TEST_CASE("SweepQuality_HDF5RoundTrip", "[HDF5][Full][Diagnostics]") {
// The per-image codes have to survive the writer and come back out of the reader. Without an
// assertion here the field can ship as all-zeros without anyone noticing.
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(6).Compression(CompressionAlgorithm::NO_COMPRESSION)
.FilePrefix("sweep_quality_roundtrip");
x.SetFileWriterFormat(FileWriterFormat::NXmxIntegrated).OverwriteExistingFiles(true);
// Images 2-3 out of beam, image 5 dead from radiation damage; the rest in no flagged range.
const std::vector<uint8_t> expected{0, 0,
static_cast<uint8_t>(SweepQualityReason::CrystalOutOfBeam) + 1,
static_cast<uint8_t>(SweepQualityReason::CrystalOutOfBeam) + 1,
0,
static_cast<uint8_t>(SweepQualityReason::RadiationDamage) + 1};
{
RegisterHDF5Filter();
StartMessage start_message;
x.FillMessage(start_message);
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
end_message.sweep_quality = expected;
end_message.sweep_quality_reasons = ReasonVocabulary();
FileWriter writer(start_message);
std::vector<int16_t> image(x.GetPixelsNum(), 42);
for (int i = 0; i < x.GetImageNum(); i++) {
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(writer.Write(message));
}
writer.WriteHDF5(end_message);
writer.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("sweep_quality_roundtrip_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->sweep_quality == expected);
// The vocabulary travels with the codes, so a consumer can name them without this source.
CHECK(dataset->sweep_quality_reasons == ReasonVocabulary());
}
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
remove("sweep_quality_roundtrip_master.h5");
}
// The command line, the wall time and the GPUs say how the result was produced, what it cost and what
// it ran on, so the report can be read on its own once the shell history is gone. They come from the
// CLI, which is the only caller that knows them; a caller that does not - the library, the viewer -
// must get a report without them rather than one claiming the run was invoked by nobody, took no
// time, and saw no GPU.
TEST_CASE("ResultReport_ProvenanceKeys", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(600);
ProcessResult result;
result.images_processed = 600;
result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
RunProvenance provenance;
provenance.command_line = "rugnux -o prefix in.h5";
provenance.wall_time_s = 262.409;
provenance.gpu_count = 4;
provenance.gpu_description = "4x NVIDIA A100-SXM4-80GB";
const auto text = RenderResultReport("prefix", "in.h5", x, result, provenance);
CHECK(text.find("\nCOMMAND_LINE= rugnux -o prefix in.h5\n") != std::string::npos);
CHECK(text.find("\nWALL_TIME= 262.41\n") != std::string::npos);
CHECK(text.find("\nGPU_COUNT= 4\n") != std::string::npos);
CHECK(text.find("\nGPU= 4x NVIDIA A100-SXM4-80GB\n") != std::string::npos);
// A machine with no GPU says so - GPU_COUNT= 0 is a statement about why the run took as long as
// it did, and only the name list has nothing to report.
RunProvenance cpu_only;
cpu_only.gpu_count = 0;
const auto cpu_text = RenderResultReport("prefix", "in.h5", x, result, cpu_only);
CHECK(cpu_text.find("\nGPU_COUNT= 0\n") != std::string::npos);
CHECK(cpu_text.find("\nGPU= ") == std::string::npos);
const auto without = RenderResultReport("prefix", "in.h5", x, result);
CHECK(without.find("COMMAND_LINE=") == std::string::npos);
CHECK(without.find("WALL_TIME=") == std::string::npos);
CHECK(without.find("GPU_COUNT=") == std::string::npos);
}
// A model that was asked for gets a section either way: the R-factors when it worked, and why it did
// not when it did not. A run that never asked for one has no section at all.
TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(600);
ProcessResult result;
result.images_processed = 600;
result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
CHECK(RenderResultReport("p", "in.h5", x, result).find("MODEL VALIDATION") == std::string::npos);
ModelValidationResult good;
good.ok = true;
good.model_path = "model.cif";
good.model_space_group_number = 96;
good.r_work = 0.1823;
good.r_free = 0.2145;
good.n_work = 16682;
good.n_free = 928;
good.mean_atom_density_sigma = 1.42;
good.fit_tested = true;
good.model_fits = true;
good.null_replicates = 5;
good.null_r_work_mean = 0.6031;
good.null_r_work_sd = 0.0094;
good.r_work_sigma = 4.12;
good.model_enantiomorph_candidate = true;
good.adopted_model_enantiomorph = true;
result.model_validation = good;
result.space_group = *gemmi::find_spacegroup_by_number(96);
const auto text = RenderResultReport("p", "in.h5", x, result);
CHECK(text.find("5. MODEL VALIDATION") != std::string::npos);
CHECK(text.find("\nMODEL_FILE= model.cif\n") != std::string::npos);
CHECK(text.find("\nR_FREE= 0.2145\n") != std::string::npos);
CHECK(text.find("\nR_FREE_REFLECTIONS= 928\n") != std::string::npos);
CHECK(text.find("\nMODEL_VALIDATION= PERFORMED\n") != std::string::npos);
CHECK(text.find("MODEL_VALIDATION= NOT_PERFORMED") == std::string::npos);
CHECK(text.find("\nMODEL_FIT= ACCEPTED\n") != std::string::npos);
// The null's internals are the gate's own workings; MODEL_FIT and MODEL_FIT_SIGMA carry the
// answer they were computed for and stay in the default report.
CHECK(text.find("MODEL_FIT_STATISTIC=") == std::string::npos);
CHECK(text.find("MODEL_FIT_NULL_REPLICATES=") == std::string::npos);
{
RunProvenance dev;
dev.developer = true;
const auto dev_text = RenderResultReport("p", "in.h5", x, result, dev);
CHECK(dev_text.find("\nMODEL_FIT_STATISTIC= R_WORK\n") != std::string::npos);
CHECK(dev_text.find("\nMODEL_FIT_NULL_REPLICATES= 5\n") != std::string::npos);
}
CHECK(text.find("\nMODEL_FIT_SIGMA= +4.12\n") != std::string::npos);
CHECK(text.find("\nMODEL_DECISIONS_TAKEN= ENANTIOMORPH\n") != std::string::npos);
// The hand is ASSUMED from the model, never determined: merged intensities cannot see it.
CHECK(text.find("\nSPACE_GROUP_ENANTIOMORPH= ASSUMED_FROM_MODEL\n") != std::string::npos);
// The same model, not accepted by the data: it still has R-factors and maps, but it decides
// nothing - the reflections stay in the group and the indexing they were merged in.
ModelValidationResult rejected = good;
rejected.model_fits = false;
rejected.r_work_sigma = 0.59;
rejected.adopted_model_enantiomorph = false;
result.model_validation = rejected;
const auto rejected_text = RenderResultReport("p", "in.h5", x, result);
CHECK(rejected_text.find("\nMODEL_FIT= REJECTED\n") != std::string::npos);
CHECK(rejected_text.find("\nMODEL_DECISIONS_TAKEN= NONE\n") != std::string::npos);
CHECK(rejected_text.find("\nR_WORK= 0.1823\n") != std::string::npos);
CHECK(rejected_text.find("\nSPACE_GROUP_ENANTIOMORPH= ASSUMED_FROM_MODEL\n") == std::string::npos);
CHECK(rejected_text.find("byte for byte") != std::string::npos);
// A model that claims nothing - already in the data's group, no merohedral ambiguity - is never
// put to the null, and NOT_TESTED has to be distinguishable from a model the data refused.
ModelValidationResult untested = good;
untested.fit_tested = false;
untested.model_fits = false;
untested.model_enantiomorph_candidate = false;
untested.adopted_model_enantiomorph = false;
result.model_validation = untested;
const auto untested_text = RenderResultReport("p", "in.h5", x, result);
CHECK(untested_text.find("\nMODEL_FIT= NOT_TESTED\n") != std::string::npos);
CHECK(untested_text.find("\nMODEL_DECISIONS_TAKEN= NONE\n") != std::string::npos);
CHECK(untested_text.find("\nR_WORK= 0.1823\n") != std::string::npos);
CHECK(untested_text.find("MODEL_FIT_NULL_MEAN=") == std::string::npos);
CHECK(untested_text.find("was tried and REJECTED") == std::string::npos);
result.model_validation = good;
ModelValidationResult failed;
failed.model_path = "broken.pdb";
failed.failure_reason = "model broken.pdb has no atoms or no unit cell";
result.model_validation = failed;
const auto failed_text = RenderResultReport("p", "in.h5", x, result);
CHECK(failed_text.find("\nMODEL_VALIDATION= NOT_PERFORMED\n") != std::string::npos);
CHECK(failed_text.find("has no atoms") != std::string::npos);
CHECK(failed_text.find("R_FREE=") == std::string::npos);
}
// get_gpu_description collapses repeats, so a four-card machine reads as one line rather than the
// same name four times. Build-independent: without CUDA there are no names and it is empty.
TEST_CASE("ResultReport_GpuDescription", "[Diagnostics]") {
const auto names = get_gpu_names();
CHECK(names.size() == static_cast<size_t>(std::max(0, get_gpu_count())));
const auto description = get_gpu_description();
if (names.empty())
CHECK(description.empty());
else
CHECK(description.find(names.front()) != std::string::npos);
}
// The verdict is the first evaluative line in the file, and it is composed from the warnings the
// sections produced - which is only possible because the report is built as a document and rendered
// afterwards. Its vocabulary is closed, so a consumer switches on it.
TEST_CASE("ResultReport_Verdict", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
// No lattice at all: the run failed, and the report says so on the first screen rather than in
// its last line.
ProcessResult nothing;
nothing.images_processed = 100;
nothing.indexing_rate = 0.0f;
const auto failed = RenderResultReport("p", "in.h5", x, nothing);
CHECK(failed.find("\nVERDICT= FAILED\n") != std::string::npos);
CHECK(failed.find("\nPATHOLOGY_FLAGS= NO_LATTICE\n") != std::string::npos);
CHECK(failed.find("WARNING: No image indexed") != std::string::npos);
// The verdict comes before the section that produced the warning behind it.
CHECK(failed.find("VERDICT=") < failed.find("LATTICE_FOUND="));
// A --mode scale run has a lattice from its input file and never fills indexing_rate. The
// "No image indexed" warning used to fire on every one of them, on reports whose own
// LATTICE_FOUND said TRUE - a warning that is false on a whole run mode is how users learn to
// skip the section.
ProcessResult rescaled;
rescaled.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
const auto scale_text = RenderResultReport("p", "in.h5", x, rescaled);
CHECK(scale_text.find("\nLATTICE_FOUND= TRUE\n") != std::string::npos);
CHECK(scale_text.find("No image indexed") == std::string::npos);
CHECK(scale_text.find("INDEXING_RATE=") == std::string::npos);
}
// Two conditions that left every other number in the report looking ordinary: a merge whose
// observations do not correlate, and a merge that measured too little of reciprocal space.
TEST_CASE("ResultReport_UnusableMergeAndCompleteness", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
ProcessResult result;
result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
result.has_merge_statistics = true;
auto &o = result.merge_statistics.overall;
o.d_max = 50.0f;
o.d_min = 2.00f;
o.unique_reflections = 1000;
o.possible_unique_reflections = 1100;
o.total_observations = 4000;
o.mean_i_over_sigma = 12.0; // strong...
o.cc_half = -0.005; // ...and uncorrelated: not equivalent reflections
o.r_meas = 1.09;
const auto text = RenderResultReport("p", "in.h5", x, result);
CHECK(text.find("\nVERDICT= UNUSABLE\n") != std::string::npos);
CHECK(text.find("WARNING: These merged data are not usable") != std::string::npos);
CHECK(text.find("UNUSABLE_MERGE") != std::string::npos);
// Rejecting more observations than were kept is the other way a merge throws the data away, and
// it improves every other statistic in the block while it does so.
ProcessResult rejected = result;
rejected.merge_statistics.overall.cc_half = 0.98;
rejected.merge_statistics.n_observations_rejected = 5000;
CHECK(RenderResultReport("p", "in.h5", x, rejected).find("threw away") != std::string::npos);
// Completeness is judged inside the fitted resolution, not over the deliberately generous cut,
// so a corner-limited detector on good data does not fire it.
ProcessResult thin = result;
thin.merge_statistics.overall.cc_half = 0.98;
thin.resolution_fit_A = 2.0;
MergeStatisticsShell shell;
shell.d_max = 50.0f;
shell.d_min = 2.0f;
shell.cc_half = 0.98;
shell.unique_reflections = 300;
shell.possible_unique_reflections = 2000;
thin.merge_statistics.shells.push_back(shell);
const auto thin_text = RenderResultReport("p", "in.h5", x, thin);
CHECK(thin_text.find("LOW_COMPLETENESS") != std::string::npos);
CHECK(thin_text.find("\nVERDICT= WARNINGS\n") != std::string::npos);
}
// FITTED_RESOLUTION is presented as "the number to quote", so it must not be quotable when the CC1/2
// curve it was read off never behaved like a fall-off.
TEST_CASE("ResultReport_FittedResolutionSuppressed", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
ProcessResult result;
result.consensus_cell = UnitCell{.a = 7.0f, .b = 7.0f, .c = 7.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
result.has_merge_statistics = true;
result.merge_statistics.overall.d_max = 6.0f;
result.merge_statistics.overall.d_min = 0.80f;
result.merge_statistics.overall.unique_reflections = 500;
result.merge_statistics.overall.possible_unique_reflections = 520;
result.merge_statistics.overall.total_observations = 2000;
result.merge_statistics.overall.cc_half = 0.90;
result.merge_statistics.overall.mean_i_over_sigma = 8.0;
// Signal to 0.90 A and nothing beyond it...
MergeStatisticsShell inner;
inner.d_max = 6.0f;
inner.d_min = 0.90f;
inner.cc_half = 0.95;
inner.unique_reflections = 250;
inner.possible_unique_reflections = 260;
MergeStatisticsShell outer = inner;
outer.d_max = 0.90f;
outer.d_min = 0.80f;
outer.cc_half = 0.02;
result.merge_statistics.shells.push_back(inner);
result.merge_statistics.shells.push_back(outer);
// ...and a fit that claims 0.70 A, finer than any shell that still reaches the target.
result.resolution_fit_A = 0.70;
const auto text = RenderResultReport("p", "in.h5", x, result);
CHECK(text.find("FITTED_RESOLUTION=") == std::string::npos);
CHECK(text.find("No FITTED_RESOLUTION is given") != std::string::npos);
CHECK(text.find("RESOLUTION_FIT") != std::string::npos);
// A fit the shells DO support is written as before.
result.resolution_fit_A = 0.95;
CHECK(RenderResultReport("p", "in.h5", x, result).find("\nFITTED_RESOLUTION= 0.95\n")
!= std::string::npos);
}
// The completeness of the range where the signal is, beside the completeness of the range that was
// written. Users compare the single overall figure against another program's and conclude the run
// lost data, when it is the deliberately generous cut diluting the denominator - so the report says
// both. But only when they are different numbers: on a detector-limited crystal every shell still
// carries signal, and printing "82.4 % to 1.22 A; 82.4 % over the full written range" says one
// thing twice.
TEST_CASE("ResultReport_CompletenessTwoRanges", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
ProcessResult result;
result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
result.has_merge_statistics = true;
auto &o = result.merge_statistics.overall;
o.d_max = 50.0f;
o.d_min = 1.20f;
o.unique_reflections = 900;
o.possible_unique_reflections = 1200;
o.total_observations = 9000;
o.cc_half = 0.95;
o.mean_i_over_sigma = 10.0;
MergeStatisticsShell inner; // signal here
inner.d_max = 50.0f;
inner.d_min = 1.50f;
inner.cc_half = 0.95;
inner.unique_reflections = 700;
inner.possible_unique_reflections = 720;
MergeStatisticsShell outer = inner; // written, but past the fall-off
outer.d_max = 1.50f;
outer.d_min = 1.20f;
outer.cc_half = 0.02;
outer.unique_reflections = 200;
outer.possible_unique_reflections = 480;
result.merge_statistics.shells.push_back(inner);
result.merge_statistics.shells.push_back(outer);
// The signal stops at 1.50 A and the reflections were written to 1.20 A: two different ranges,
// so both numbers are reported and the outer shells are visibly what dilutes the overall one.
const auto narrower = RenderResultReport("p", "in.h5", x, result);
CHECK(narrower.find("% to 1.50 A;") != std::string::npos);
CHECK(narrower.find("\nCOMPLETENESS= 75.0\n") != std::string::npos);
// Detector-limited: every shell still reaches the target, so the range where the signal is IS
// the written range and the second number would repeat the first.
result.merge_statistics.shells[1].cc_half = 0.95;
result.merge_statistics.shells[1].possible_unique_reflections = 200;
const auto same = RenderResultReport("p", "in.h5", x, result);
CHECK(same.find(" % to ") == std::string::npos);
CHECK(same.find("over the full written range") != std::string::npos);
}
// The TNCS_* block is an interface in the same way the keys above are, and the distinction between
// FALSE, INCONCLUSIVE and NOT_MEASURED is the whole reason it is an enum rather than a boolean: a
// merge the detector could not run on must not read as a crystal that was tested and found clean.
TEST_CASE("ResultReport_TranslationalNCS", "[Diagnostics]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(600);
ProcessResult result;
result.images_processed = 600;
result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
result.space_group = *gemmi::find_spacegroup_by_number(96);
result.has_merge_statistics = true;
{
const auto text = RenderResultReport("prefix", "in.h5", x, result);
CHECK(text.find("\nTNCS_DETECTED= NOT_MEASURED\n") != std::string::npos);
CHECK(text.find("not a\n statement that this crystal has none") != std::string::npos);
}
result.tncs.measurable = true;
result.tncs.n_reflections = 4000;
result.tncs.peak_percent = 54.2;
result.tncs.peak_z = 42.3;
result.tncs.null_mean = 9.9;
result.tncs.null_sd = 1.05;
result.tncs.modulation_measured = true;
result.tncs.vector_frac = {0.0, 0.5, 0.25};
result.tncs.vector_length_A = 40.7;
result.tncs.modulation = 12.50;
result.tncs.modulation_null = 1.48;
result.tncs.detected = true;
result.tncs.near_extinct_class = true;
const auto text = RenderResultReport("prefix", "in.h5", x, result);
CHECK(text.find("\nTNCS_DETECTED= TRUE\n") != std::string::npos);
CHECK(text.find("\nTNCS_PATTERSON_PEAK_PCT= 54.2\n") != std::string::npos);
CHECK(text.find("\nTNCS_PATTERSON_PEAK_Z= 42.3\n") != std::string::npos);
CHECK(text.find("\nTNCS_PATTERSON_NULL_PCT= 9.9 +- 1.05\n") != std::string::npos);
CHECK(text.find("\nTNCS_VECTOR= 0.000 0.500 0.250\n") != std::string::npos);
CHECK(text.find("\nTNCS_VECTOR_LENGTH= 40.7\n") != std::string::npos);
CHECK(text.find("\nTNCS_MODULATION= 12.50\n") != std::string::npos);
CHECK(text.find("\nTNCS_MODULATION_NULL= 1.48\n") != std::string::npos);
CHECK(text.find("\nTNCS_SUBLATTICE= NEAR_EXTINCT_CLASS\n") != std::string::npos);
// A detection is a condition a consumer switches on, not only a sentence to read.
CHECK(text.find("\nPATHOLOGY_FLAGS= PSEUDO_TRANSLATION LATTICE_TRANSLATION\n")
!= std::string::npos);
CHECK(text.find("\nVERDICT= WARNINGS\n") != std::string::npos);
// TNCS_REFLECTIONS is a pipeline internal and stays out of the default report.
CHECK(text.find("\nTNCS_REFLECTIONS=") == std::string::npos);
RunProvenance dev;
dev.developer = true;
CHECK(RenderResultReport("prefix", "in.h5", x, result, dev).find("\nTNCS_REFLECTIONS= 4000\n")
!= std::string::npos);
// A translation the data are EXACTLY invariant under is a lattice vector, and is reported as one
// rather than as a pseudo-symmetry - this is the centred-lattice-merged-in-P1 case.
result.tncs.detected = false;
result.tncs.near_extinct_class = false;
result.tncs.undeclared_lattice_translations.push_back({0.5, 0.5, 0.0});
const auto lat = RenderResultReport("prefix", "in.h5", x, result);
CHECK(lat.find("\nTNCS_DETECTED= FALSE\n") != std::string::npos);
CHECK(lat.find("\nUNDECLARED_LATTICE_TRANSLATION= 0.500 0.500 0.000\n") != std::string::npos);
CHECK(lat.find("\nPATHOLOGY_FLAGS= LATTICE_TRANSLATION\n") != std::string::npos);
CHECK(lat.find("\nTNCS_SUBLATTICE=") == std::string::npos);
}
Reference in New Issue View Git Blame Copy Permalink