mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 78 Wiki Activity
Files
rc169
Jungfraujoch/tests/ResultReportTest.cpp
T
leonarski_fandClaude Opus 5 d00175c8b4 resolution: a cut the shell table refutes is not quoted 
One run was written to 1.089 A where I/sigma was zero and R_meas 1776 per cent,
and the summary said only "Merged to 1.09 A". Two things let that through.

The logistic's crossing was taken wherever the fitted curve met the target, even
when that point lay past every bin the curve was fitted through - an extrapolation
of a fall-off the data never showed, quoted as a measurement. The crossing now has
to lie inside the fitted bins, with the one-shell extension still there to spare;
outside them the number is read off the bins themselves instead.

And the shipped guard asked only for the finest shell whose CC1/2 still reached the
target, with no requirement that the curve get there monotonically. A noise shell
that climbs back over the bar therefore became the edge of the data. It is the
climbing back that disqualifies it, and the program already noticed - it printed
"CC1/2 is not monotone with resolution" and then cut there anyway, because the test
was developer-only and decided nothing. It decides now, in the report everyone
reads, and the duplicate is gone so there is one such test rather than two.

The dataset above is written to 1.527 A: CC1/2 0.209 to 0.688, I/sigma 0.40 to
1.26, R_meas 245 to 145 per cent. Forty-four of fifty-one reports are unchanged to
the character; of the seven that move, five are cut coarser and every headline
number of all five improves, one loses a fit that was an extrapolation without
changing anything written, and one gains a quotable fit and loses a warning.

This is a guard, not the cause. On four of those five the reflections doing the
damage are ice, which the resolution fit now leaves out for its own reasons; the
guard still earns its place, because on those four removing the ice alone does not
stop the cut being quoted past what the shells support.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011GxZqDiFP3KqriBhNdcR56
2026-09-14 14:17:39 +02:00

710 lines
37 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,
.disposition = FrameDisposition::Rejected,
.delta_cc_half = -0.0312f, .delta_cc_half_se = 0.0060f});
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,
.disposition = FrameDisposition::Downgraded,
.delta_cc_half = 0.0011f, .delta_cc_half_se = 0.0042f});
out.frames_merged = 450;
out.frames_downgraded = 100;
out.frames_rejected = 50;
out.rejected_deg = 5.0f;
return out;
}
std::vector<std::string> ReasonVocabulary() {
std::vector<std::string> out;
for (int r = 0; r <= static_cast<int>(SweepQualityReason::InconsistentWithMerge); ++r)
out.emplace_back(SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
return out;
}
std::vector<std::string> DispositionVocabulary() {
std::vector<std::string> out;
for (int d = 0; d <= static_cast<int>(FrameDisposition::Rejected); ++d)
out.emplace_back(FrameDispositionCode(static_cast<FrameDisposition>(d)));
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() == 6);
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");
CHECK(codes[5] == "inconsistent_with_merge");
// The disposition vocabulary is an interface for the same reasons, and the per-image HDF5 array
// indexes it from 0.
CHECK(std::string(FrameDispositionCode(FrameDisposition::Merged)) == "merged");
CHECK(std::string(FrameDispositionCode(FrameDisposition::Downgraded)) == "downgraded");
CHECK(std::string(FrameDispositionCode(FrameDisposition::Rejected)) == "rejected");
}
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= 8\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 inconsistent_with_merge\n") != std::string::npos);
CHECK(dev_text.find("\nSWEEP_DISPOSITIONS= merged downgraded rejected\n") != std::string::npos);
// The disposition headline: frames AND degrees, so that re-slicing cannot move the percentage
// without the degrees moving with it.
CHECK(text.find("\nFRAMES_MERGED= 450\n") != std::string::npos);
CHECK(text.find("\nFRAMES_DOWNGRADED= 100\n") != std::string::npos);
CHECK(text.find("\nFRAMES_REJECTED= 50\n") != std::string::npos);
CHECK(text.find("\nFRAMES_REJECTED_PCT= 8.33\n") != std::string::npos);
CHECK(text.find("\nROTATION_REJECTED_DEG= 5.0\n") != std::string::npos);
CHECK(text.find("600 frames / 60.0 deg; 450 merged (75.0%), 100 downgraded (16.7%)")
!= 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 its disposition beside what keeping it costs the merged intensities.
CHECK(text.find("rejected -0.0312 0.0060") != std::string::npos);
CHECK(text.find("downgraded +0.0011 0.0042") != 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("\n# WARNING: Frames 100-149 out of beam (5.0 deg,") != std::string::npos);
CHECK(text.find("\n# WARNING: Frames 400-499 radiation damage (10.0 deg,") != std::string::npos);
CHECK(text.find("REJECTED: keeping them lowered CC1/2 by 0.0312 +/- 0.0060 over the "
"reflections they touch") != std::string::npos);
CHECK(text.find("kept at their own reduced weight; removing them would change CC1/2 by -0.0011")
!= 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};
// What was DONE about them: the out-of-beam pair is dropped, the damaged tail stays in at its own
// reduced weight, the rest merges. Indexed from 0, unlike the reason codes above.
const std::vector<uint8_t> expected_disposition{
static_cast<uint8_t>(FrameDisposition::Merged),
static_cast<uint8_t>(FrameDisposition::Merged),
static_cast<uint8_t>(FrameDisposition::Rejected),
static_cast<uint8_t>(FrameDisposition::Rejected),
static_cast<uint8_t>(FrameDisposition::Merged),
static_cast<uint8_t>(FrameDisposition::Downgraded)};
{
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();
end_message.frame_disposition = expected_disposition;
end_message.frame_disposition_codes = DispositionVocabulary();
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());
CHECK(dataset->frame_disposition == expected_disposition);
CHECK(dataset->frame_disposition_codes == DispositionVocabulary());
}
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("The table reads 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);
// CC1/2 that drops through the target and climbs back above it is not a fall-off either, and the
// finest shell reaching the target is then noise correlating with noise - so a fit the test above
// accepts must still be suppressed, and the user told so in the default report.
ProcessResult climbs = result;
climbs.merge_statistics.shells.back().cc_half = 0.05;
MergeStatisticsShell noise = climbs.merge_statistics.shells.back();
noise.d_max = 0.80f;
noise.d_min = 0.70f;
noise.cc_half = 0.60;
climbs.merge_statistics.shells.push_back(noise);
const auto climbs_text = RenderResultReport("p", "in.h5", x, climbs);
CHECK(climbs_text.find("FITTED_RESOLUTION=") == std::string::npos);
CHECK(climbs_text.find("climbs back above it at 0.70 A") != std::string::npos);
CHECK(climbs_text.find("RESOLUTION_FIT") != 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