mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 3 Packages Projects Releases 74 Wiki Activity
Files
rc167
Jungfraujoch/tests/ResultReportTest.cpp
T
leonarski_fandClaude Opus 5 004fa5a781 model validation: the model says, shell by shell, where the data still carry signal 
A placed and scaled model gives an expected intensity for every reflection, and comparing that
with the merged intensity is a measurement nothing else in the report can make: CC1/2 and R_meas
compare the data with themselves, and can only say whether a shell is reproducible, not whether
it is right. CC(model, data) is reported per shell, on the merge's own shells so the rows line
up, with the reflection count and a significance beside each.

Read it in one direction only. A correlation significantly above zero proves signal is present,
because a model cannot invent agreement with data it never saw. A correlation near zero proves
nothing, since the model may be the thing at fault - one crystal in this corpus reads negative
across its outer shells against a deposited model of a different crystal form. So this can
support keeping a shell and must never be used to discard one, and nothing in the pipeline acts
on it.

Intensities rather than amplitudes, matching CC1/2 and CCref beside it. The raw merged intensity
rather than the French-Wilson amplitude: that amplitude is a posterior mean under a Wilson prior
which pulls weak reflections towards their shell mean, and the shells this number is read in are
exactly the weak ones, so it would correlate the prior with itself.

The free flag is ignored on purpose. Nothing has been refined against these reflections, so the
correlation is unbiased over all of a shell rather than over the test set alone - which is where
the outer-shell reflection counts are, and with them the difference between a statistic that can
resolve something and one that cannot.

REPORT_VERSION is 7.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
2026-09-07 13:12:32 +02:00

642 lines
33 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 <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= 7\n") != std::string::npos);
// SOHNCKE_SPACE_GROUP names the best group a chiral crystal could have, and it comes from the
// space-group SEARCH. This fixture is given its group rather than searching for one, so there is
// no Sohncke candidate to name and the key is absent - which is the honest behaviour and the
// reason the changelog says "every run that determined a space group" and not "every run".
CHECK(text.find("\nSOHNCKE_SPACE_GROUP=") == 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);
// A run that measured no CC(model, data) writes no key for it.
CHECK(text.find("CC_MODEL_OVERALL=") == std::string::npos);
CHECK(text.find("CC_MODEL_CONFIRMED_TO_D_MIN=") == 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);
// CC(model, data) by shell, with the deepest shell short of significance: the confirmed limit is
// the last shell that reached it, and it is a lower bound - the shell past it is not thereby empty.
ModelValidationResult with_cc = good;
with_cc.cc_model_shells = {{3.20f, 0.9412, 4210, 91.4}, {2.10f, 0.5533, 3980, 40.2},
{1.80f, 0.0412, 2110, 1.9}};
with_cc.cc_model_overall = 0.8123;
with_cc.cc_model_n = 10300;
result.model_validation = with_cc;
const auto cc_text = RenderResultReport("p", "in.h5", x, result);
CHECK(cc_text.find("\nCC_MODEL_OVERALL= 0.8123\n") != std::string::npos);
CHECK(cc_text.find("\nCC_MODEL_REFLECTIONS= 10300\n") != std::string::npos);
CHECK(cc_text.find("\nCC_MODEL_CONFIRMED_TO_D_MIN= 2.10\n") != std::string::npos);
CHECK(cc_text.find("D_MIN CC_MODEL N SIGMA") != std::string::npos);
CHECK(cc_text.find(" 1.80 0.0412 2110 +1.9") != std::string::npos);
CHECK(cc_text.find("Read it in ONE direction only") != std::string::npos);
// Nothing significant anywhere: the key still has to be written, saying so.
ModelValidationResult no_signal = with_cc;
for (auto &sh : no_signal.cc_model_shells)
sh.sigma = 1.0;
result.model_validation = no_signal;
CHECK(RenderResultReport("p", "in.h5", x, result)
.find("\nCC_MODEL_CONFIRMED_TO_D_MIN= NONE\n") != 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