leonarski_fandClaude Opus 5
54adcaafcc
The per-image resolution estimate was the 5th percentile of the spot d spacings - an extreme order statistic, so it measured where detection stops rather than how well the crystal diffracts. A large cell puts more reflections past the same threshold and scored better than a small cell that diffracts further; intensity was not used at all, so a weak crystal padded with spurious high-resolution detections ran away; and nothing clamped the answer to what the detector can deliver. Against the resolution the merged data actually reach it was 42% out in log-RMS, with 1 of 38 rotation datasets inside 0.2 A. Take instead the 1/d^2 beyond which 30% of the sum of sqrt(I) over the image's non-ice spots lies, report 1/(2.25 sqrt of it), clamp at the detector corner, and take the median over images. A quantile from the middle of the distribution measures the shape of the falloff - the crystal's own exp(-B/2d^2) - where an extreme one measures the threshold. sqrt(I) is the Poisson significance of a summed photon count, so a marginal high-resolution detection cannot carry the answer and neither can a handful of very strong low-resolution reflections. The 2.25 is the multiplicity gain: merging keeps measuring intensities a fixed factor in 1/d past the point where a single frame detects them. Spearman 0.881 -> 0.954, log-RMS 42% -> 8.9%, median error 0.79 -> 0.07 A, and 32 of 38 within 0.2 A. Both constants sit on a broad plateau, the scale is stable across dataset halves and across resolution ranges, and no second predictor survives leave-one-out. The residual is around 9%, set by multiplicity, symmetry and radiation damage - none of which a spot list can see. The estimate feeds only reporting: the image stream, HDF5, the plots, the scan result and the preview ring. It sets no cutoff and no search limit, and the scaling and merging output is byte-identical. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016NNnL26LAvruQ9eLUUWvrJ
303 lines
17 KiB
C++
303 lines
17 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <fstream>
|
|
#include <sstream>
|
|
|
|
#include <spdlog/fmt/fmt.h>
|
|
|
|
#include "../common/GitInfo.h"
|
|
#include "../common/time_utc.h"
|
|
#include "../image_analysis/scale_merge/Merge.h"
|
|
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
|
|
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
|
|
|
|
#include "ResultReport.h"
|
|
|
|
namespace {
|
|
// The version of this file format. Bumped when a key is renamed or removed, a table column moves,
|
|
// or a reason code changes meaning - a consumer can gate on it.
|
|
constexpr int REPORT_VERSION = 1;
|
|
|
|
const char *BANNER = " ******************************************************************************";
|
|
|
|
void Section(std::ostream &os, const std::string &title) {
|
|
os << "\n" << BANNER << "\n " << title << "\n" << BANNER << "\n\n";
|
|
}
|
|
|
|
// Every number a consumer might want is written as one of these, so it is one grep away.
|
|
template <class T> void Key(std::ostream &os, const char *key, const T &value) {
|
|
os << key << "= " << value << "\n";
|
|
}
|
|
|
|
std::string CellString(const UnitCell &c) {
|
|
return fmt::format("{:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}", c.a, c.b, c.c,
|
|
c.alpha, c.beta, c.gamma);
|
|
}
|
|
}
|
|
|
|
std::string RenderResultReport(const std::string &output_prefix,
|
|
const std::string &input_file,
|
|
const DiffractionExperiment &experiment,
|
|
const ProcessResult &result) {
|
|
std::ostringstream os;
|
|
const bool rotation = experiment.IsRotationIndexing();
|
|
const bool merged = result.has_merge_statistics;
|
|
std::vector<std::string> warnings = result.warnings;
|
|
|
|
os << BANNER << "\n"
|
|
<< " RUGNUX PROCESSING REPORT\n"
|
|
<< BANNER << "\n\n"
|
|
<< " What this run determined, written next to its other output. The `KEY= value` lines and\n"
|
|
<< " the tables below are a stable interface - a script greps them, and REPORT_VERSION says\n"
|
|
<< " when that interface last changed. Timing, rates and per-image progress are not here;\n"
|
|
<< " they are on stdout.\n\n";
|
|
|
|
Key(os, "REPORT_VERSION", REPORT_VERSION);
|
|
Key(os, "RUGNUX_VERSION", jfjoch_version());
|
|
if (!jfjoch_git_sha1().empty())
|
|
Key(os, "RUGNUX_GIT", jfjoch_git_sha1().substr(0, 6) + " " + jfjoch_git_date());
|
|
Key(os, "DATE", time_UTC(std::chrono::system_clock::now()));
|
|
Key(os, "INPUT_FILE", input_file);
|
|
Key(os, "OUTPUT_PREFIX", output_prefix);
|
|
|
|
// ---------------------------------------------------------------- 1. DATA SET
|
|
Section(os, "1. DATA SET");
|
|
Key(os, "EXPERIMENT_TYPE", rotation ? "ROTATION" : "STILLS");
|
|
Key(os, "IMAGES_PROCESSED", result.images_processed);
|
|
Key(os, "WAVELENGTH", fmt::format("{:.5f}", experiment.GetWavelength_A()));
|
|
if (const auto gonio = experiment.GetGoniometer()) {
|
|
Key(os, "OSCILLATION_RANGE", fmt::format("{:.4f}", gonio->GetIncrement_deg()));
|
|
Key(os, "STARTING_ANGLE", fmt::format("{:.3f}", gonio->GetStart_deg()));
|
|
const auto ax = gonio->GetAxis();
|
|
Key(os, "ROTATION_AXIS", fmt::format("{:.6f} {:.6f} {:.6f}", ax.x, ax.y, ax.z));
|
|
}
|
|
Key(os, "DETECTOR_DISTANCE", fmt::format("{:.3f}", result.used_distance_mm));
|
|
Key(os, "BEAM_CENTRE", fmt::format("{:.2f} {:.2f}", result.used_beam_x_pxl, result.used_beam_y_pxl));
|
|
os << "\n"
|
|
<< " The distance and beam centre above are the ones this result was integrated at, which on\n"
|
|
<< " a rotation run is the post-refined geometry rather than the values in the input file.\n";
|
|
|
|
if (result.spot_resolution_estimate_A.has_value()) {
|
|
os << "\n";
|
|
Key(os, "SPOT_RESOLUTION_ESTIMATE", fmt::format("{:.2f}", *result.spot_resolution_estimate_A));
|
|
os << "\n"
|
|
<< " How far the merged data are expected to reach, read off the found spots alone - no\n"
|
|
<< " lattice, no integration, no merge. It is a prediction, good to about 0.2 A on the\n"
|
|
<< " rotation data it was calibrated on, and it is not what the run achieved: compare it\n"
|
|
<< " with INCLUDE_RESOLUTION_RANGE in section 5.\n";
|
|
}
|
|
|
|
if (result.pass_count > 1) {
|
|
os << "\n";
|
|
Key(os, "PASS", fmt::format("{} of {}", result.pass_number, result.pass_count));
|
|
Key(os, "PASS_DECISION", result.pass_decision);
|
|
os << "\n"
|
|
<< " A rotation run integrates twice: once at the geometry in the input file, then again at\n"
|
|
<< " the post-refined geometry. Every number in this report describes the pass named above,\n"
|
|
<< " whose files are " << output_prefix << ".*; the header-geometry pass is not written.\n";
|
|
}
|
|
|
|
// ---------------------------------------------------------------- 2. INDEXING
|
|
Section(os, "2. INDEXING");
|
|
if (result.indexing_rate.has_value())
|
|
Key(os, "INDEXING_RATE", fmt::format("{:.4f}", result.indexing_rate.value()));
|
|
Key(os, "LATTICE_FOUND", (result.consensus_cell.has_value() ? "TRUE" : "FALSE"));
|
|
if (result.consensus_cell.has_value())
|
|
Key(os, "UNIT_CELL_CONSTANTS", CellString(*result.consensus_cell));
|
|
if (result.space_group_number.has_value())
|
|
Key(os, "SPACE_GROUP_NUMBER", result.space_group_number.value());
|
|
if (result.indexing_rate.value_or(0.0f) <= 0.0f)
|
|
warnings.emplace_back("No image indexed - no crystal lattice was determined from this dataset");
|
|
|
|
// ---------------------------------------------- 3. GEOMETRY POST-REFINEMENT
|
|
if (result.post_refine.has_value()) {
|
|
const auto &pr = *result.post_refine;
|
|
Section(os, "3. GEOMETRY POST-REFINEMENT");
|
|
os << " The rotation two-pass fits the detector distance and beam centre from the observed spot\n"
|
|
<< " positions, and the cell scale and rotation axis from the observed rocking angles. Each\n"
|
|
<< " step is committed only if it improves a held-out residual.\n\n";
|
|
Key(os, "POSTREFINE_EVENTS_USED", pr.events_used);
|
|
Key(os, "POSTREFINE_OBS_USED", pr.obs_used);
|
|
Key(os, "POSTREFINE_CELL_COMMITTED", pr.cell_refined ? "TRUE" : "FALSE");
|
|
Key(os, "POSTREFINE_DETECTOR_COMMITTED", pr.detector_refined ? "TRUE" : "FALSE");
|
|
Key(os, "POSTREFINE_DISTANCE", fmt::format("{:.3f} -> {:.3f}", pr.distance_before_mm,
|
|
pr.distance_after_mm));
|
|
Key(os, "POSTREFINE_BEAM_CENTRE", fmt::format("{:.2f} {:.2f} -> {:.2f} {:.2f}",
|
|
pr.beam_x_before_px, pr.beam_y_before_px,
|
|
pr.beam_x_after_px, pr.beam_y_after_px));
|
|
Key(os, "GONIOMETER_ROTATION_SCALE", fmt::format("{:.5f}", pr.rotation_scale));
|
|
Key(os, "GONIOMETER_ROTATION_SCALE_SUSPECT", pr.rotation_scale_suspect ? "TRUE" : "FALSE");
|
|
os << "\n GONIOMETER_ROTATION_SCALE is the factor by which the stage actually turned relative to\n"
|
|
<< " the angles stored in the file (which are the commanded ones). 1.0 = they agree. It drives\n"
|
|
<< " the second integration pass only when SUSPECT is TRUE - both cross-validated and outside\n"
|
|
<< " the tolerance - since a stage that is in fact well calibrated must be left alone. A\n"
|
|
<< " manual --rotation-scale replaces it and is applied to both passes.\n";
|
|
if (pr.rotation_scale_suspect)
|
|
warnings.emplace_back(fmt::format(
|
|
"The goniometer turned by a factor {:.5f} of the angles stored in the file - the "
|
|
"stage rotation looks mis-calibrated by {:+.2f}%. The correction was applied to this "
|
|
"run, but the fault is in the hardware and should be fixed there",
|
|
pr.rotation_scale, 100.0 * (pr.rotation_scale - 1.0)));
|
|
}
|
|
|
|
// ---------------------------------------------- 4. SPACE GROUP DETERMINATION
|
|
Section(os, "4. SPACE GROUP DETERMINATION");
|
|
if (result.space_group_search.has_value()) {
|
|
Key(os, "SPACE_GROUP_SEARCH", "DE_NOVO");
|
|
os << "\n" << SearchSpaceGroupResultToText(*result.space_group_search) << "\n";
|
|
} else if (result.space_group_number.has_value()) {
|
|
Key(os, "SPACE_GROUP_SEARCH", "FIXED");
|
|
os << "\n The space group was given, not determined here.\n";
|
|
} else {
|
|
Key(os, "SPACE_GROUP_SEARCH", "NONE");
|
|
os << "\n No space group was determined.\n";
|
|
}
|
|
|
|
// ---------------------------------------------------- 5. SCALING AND MERGING
|
|
Section(os, "5. SCALING AND MERGING");
|
|
if (!merged) {
|
|
Key(os, "MERGE", "NOT_PERFORMED");
|
|
os << "\n No scaling or merging was performed on this run, so there are no merging statistics, no\n"
|
|
<< " error model, and no sweep-quality diagnosis below. The integrated reflections are in\n"
|
|
<< " " << output_prefix << "_process.h5.\n";
|
|
} else {
|
|
const auto &o = result.merge_statistics.overall;
|
|
Key(os, "MERGE", "PERFORMED");
|
|
Key(os, "INCLUDE_RESOLUTION_RANGE", fmt::format("{:.3f} {:.3f}", o.d_max, o.d_min));
|
|
Key(os, "FRIEDELS_LAW", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE");
|
|
Key(os, "UNIQUE_REFLECTIONS", o.unique_reflections);
|
|
Key(os, "TOTAL_OBSERVATIONS", o.total_observations);
|
|
Key(os, "COMPLETENESS", o.possible_unique_reflections > 0
|
|
? fmt::format("{:.1f}", 100.0 * o.unique_reflections / o.possible_unique_reflections)
|
|
: std::string("nan"));
|
|
Key(os, "MULTIPLICITY", o.unique_reflections > 0
|
|
? fmt::format("{:.2f}", static_cast<double>(o.total_observations) / o.unique_reflections)
|
|
: std::string("nan"));
|
|
Key(os, "I_OVER_SIGMA", fmt::format("{:.2f}", o.mean_i_over_sigma));
|
|
Key(os, "R_MEAS", fmt::format("{:.4f}", o.r_meas));
|
|
Key(os, "CC_HALF", fmt::format("{:.4f}", o.cc_half));
|
|
Key(os, "SIGANO", fmt::format("{:.3f}", o.abs_diff_over_sigma_anomalous));
|
|
Key(os, "WILSON_B", fmt::format("{:.2f}", result.merge_statistics.wilson_b));
|
|
// The error model in XDS's convention, so the numbers are directly comparable with a CORRECT.LP.
|
|
Key(os, "ERROR_MODEL_A", fmt::format("{:.4f}", result.error_model_a));
|
|
Key(os, "ERROR_MODEL_B", fmt::format("{:.4e}", result.error_model_b));
|
|
Key(os, "ISA", fmt::format("{:.2f}", result.error_model_isa));
|
|
if (result.error_model_isa_asymptotic > 0.0)
|
|
Key(os, "ISA_ASYMPTOTIC", fmt::format("{:.2f}", result.error_model_isa_asymptotic));
|
|
Key(os, "REFERENCE_DATA_USED", result.has_reference ? "TRUE" : "FALSE");
|
|
// The shell table straight off the statistics rather than result.merge_statistics_text: that
|
|
// string also carries the twinning analysis and the advisories, which have sections of their own.
|
|
os << "\n ERROR_MODEL_A / ERROR_MODEL_B are in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2),\n"
|
|
<< " so ISA = 1/sqrt(a*b) means what CORRECT.LP's ISa means. ISA_ASYMPTOTIC, where present,\n"
|
|
<< " is the strong-reflection tier only.\n\n"
|
|
<< result.merge_statistics;
|
|
}
|
|
|
|
// --------------------------------------------------------------- 6. TWINNING
|
|
if (merged && result.twinning.l_test_pairs > 0) {
|
|
Section(os, "6. TWINNING");
|
|
Key(os, "TWINNING_SUSPECTED", result.twinning.twinning_suspected ? "TRUE" : "FALSE");
|
|
Key(os, "L_TEST_MEAN_ABS_L", fmt::format("{:.4f}", result.twinning.mean_abs_l));
|
|
Key(os, "L_TEST_MEAN_L_SQUARED", fmt::format("{:.4f}", result.twinning.mean_l_squared));
|
|
Key(os, "SECOND_MOMENT_I", fmt::format("{:.4f}", result.twinning.second_moment));
|
|
Key(os, "ESTIMATED_TWIN_FRACTION", fmt::format("{:.3f}", result.twinning.estimated_twin_fraction));
|
|
os << "\n" << TwinningAnalysisToText(result.twinning) << "\n";
|
|
if (result.twinning.twinning_suspected)
|
|
warnings.emplace_back(fmt::format(
|
|
"Twinning is indicated (<|L|> = {:.3f}, <I^2>/<I>^2 = {:.3f}, estimated twin "
|
|
"fraction {:.2f}) - refine against the merged data with care",
|
|
result.twinning.mean_abs_l, result.twinning.second_moment,
|
|
result.twinning.estimated_twin_fraction));
|
|
}
|
|
|
|
// ------------------------------------------------------- 7. RADIATION DAMAGE
|
|
if (!result.radiation_damage_text.empty()) {
|
|
Section(os, "7. RADIATION DAMAGE");
|
|
// A number, or a word saying why there is none: NOT_A_TREND where the per-batch curve was measured
|
|
// but no straight line describes it (damage is progressive, so that curve is not dose), NOT_MEASURED
|
|
// where the monitor could not run at all.
|
|
const double db = result.merge_statistics.radiation_damage_delta_b;
|
|
Key(os, "RADIATION_DAMAGE_RELATIVE_B",
|
|
std::isfinite(db) ? fmt::format("{:.2f}", db)
|
|
: result.merge_statistics.radiation_damage_b_batch.empty() ? std::string("NOT_MEASURED")
|
|
: std::string("NOT_A_TREND"));
|
|
os << "\n" << result.radiation_damage_text << "\n";
|
|
}
|
|
|
|
// ------------------------------------------------------ 8. SWEEP QUALITY
|
|
const auto &sq = result.merge_statistics.sweep_quality;
|
|
Section(os, "8. SWEEP QUALITY");
|
|
os << " Stretches of the sweep over which the crystal delivered much less than the rest of the run.\n"
|
|
<< " REASON comes from a closed vocabulary, listed below so a consumer can tell an unknown code\n"
|
|
<< " from a missing one. SEVERITY is the fraction of the run's typical diffracting power missing\n"
|
|
<< " over the range (0 = as good as the run, 1 = nothing at all); SCALE and CC are the range's\n"
|
|
<< " mean per-image scale and CC-to-merge relative to the run median; INDEXED is the fraction of\n"
|
|
<< " the range's frames that were scaled at all. Nothing is excluded on the strength of this.\n\n";
|
|
Key(os, "SWEEP_QUALITY_STATUS", sq.measured ? "COMPUTED" : "NOT_COMPUTED");
|
|
Key(os, "SWEEP_QUALITY_COUNT", sq.ranges.size());
|
|
{
|
|
std::string codes;
|
|
for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
|
|
codes += (codes.empty() ? "" : " ")
|
|
+ std::string(SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
|
|
Key(os, "SWEEP_QUALITY_REASONS", codes);
|
|
}
|
|
if (sq.measured) {
|
|
Key(os, "SWEEP_ROTATION", fmt::format("{:.1f}", sq.sweep_deg));
|
|
Key(os, "FLUX_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.flux_peak_to_trough));
|
|
Key(os, "SCALE_MODULATION_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.modulation_peak_to_trough));
|
|
}
|
|
os << "\n"
|
|
<< " FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED\n"
|
|
<< " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n";
|
|
for (const auto &r : sq.ranges) {
|
|
os << fmt::format(" {:11d} {:11d} {:9d} {:8.1f} {:<20} {:8.2f} {:6.2f} {:6.2f} {:8.2f}\n",
|
|
r.first_image, r.last_image, r.last_image - r.first_image + 1, r.rotation_deg,
|
|
SweepQualityReasonCode(r.reason), r.severity, r.mean_relative_scale,
|
|
r.mean_relative_cc, r.indexed_fraction);
|
|
warnings.push_back(fmt::format(
|
|
"Frames {}-{} {} ({:.1f} deg, scale {:.2f} and CC {:.2f} of the run, {:.0f}% scaled)",
|
|
r.first_image, r.last_image, SweepQualityReasonText(r.reason), r.rotation_deg,
|
|
r.mean_relative_scale, r.mean_relative_cc, 100.0 * r.indexed_fraction));
|
|
}
|
|
os << " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n";
|
|
|
|
// --------------------------------------------------------------- 9. WARNINGS
|
|
if (result.cancelled)
|
|
warnings.emplace_back(fmt::format("Processing was cancelled after {} images - this report "
|
|
"describes an incomplete run", result.images_processed));
|
|
Section(os, "9. WARNINGS");
|
|
os << " Everything that needs a person's attention, one line each, marked so a script can find\n"
|
|
<< " them with a single grep for \"WARNING:\".\n\n";
|
|
Key(os, "WARNING_COUNT", warnings.size());
|
|
os << "\n";
|
|
for (const auto &w : warnings)
|
|
os << "WARNING: " << w << "\n";
|
|
if (warnings.empty())
|
|
os << " (none)\n";
|
|
|
|
os << "\n" << BANNER << "\n END OF REPORT\n" << BANNER << "\n";
|
|
return os.str();
|
|
}
|
|
|
|
void WriteResultReport(const std::string &output_prefix,
|
|
const std::string &input_file,
|
|
const DiffractionExperiment &experiment,
|
|
const ProcessResult &result,
|
|
Logger &logger) {
|
|
if (output_prefix.empty())
|
|
return; // "compute the statistics, persist nothing"
|
|
|
|
const std::string filename = output_prefix + "_report.txt";
|
|
// The report is unconditional, so it must never be the reason a run fails: a run that produced a
|
|
// good .mtz must survive an unwritable path or a full disk. Report the failure and carry on.
|
|
try {
|
|
std::ofstream file(filename);
|
|
file.exceptions(std::ios::failbit | std::ios::badbit);
|
|
file << RenderResultReport(output_prefix, input_file, experiment, result);
|
|
} catch (const std::exception &e) {
|
|
logger.Warning("Could not write the results report {}: {}", filename, e.what());
|
|
}
|
|
}
|