mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 7 Packages Projects Releases 71 Wiki Activity
Files
c8b05295c6c6b70fe52a3a52bf289a18ebb217d9
Jungfraujoch/image_analysis/WriteReflections.cpp
T
leonarski_fandClaude Opus 5 abb61e2708 Unmerged export: carry the sort key with the part here too 
SumRockingEvents assembles rocking events for the unmerged MTZ through
the same indirect comparator the anisotropy diagnostic had, on a larger
array - it has neither the scaled-image test nor the usable() filter in
front of it. It is latent, since --export-unmerged is off by default, but
it is the same defect and it is a mechanical fix.

Same treatment: the key travels with the part, plus an exact reserve.
p_unmerged.mtz, p_unmerged_partials.mtz and p.hkl are md5-identical
across eight runs on two crystals. With the export on it is worth about
0.9 s on a seventeen-million-partial run.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
2026-08-26 14:11:56 +02:00

894 lines
48 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-only
#include "WriteReflections.h"
#include "scale_merge/Merge.h"
#include "scale_merge/HKLKey.h"
#include "scale_merge/TwinningAnalysis.h"
#include "bragg_integration/SystematicAbsence.h"
#include "../common/ParallelFor.h"
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <map>
#include <set>
#include <tuple>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <stdexcept>
#include <ctime>
#include <chrono>
#include <gemmi/mtz.hpp>
#include "../common/GitInfo.h"
namespace {
/// Current date in ISO-8601 (YYYY-MM-DD) for the _audit block.
std::string CurrentDateISO() {
auto now = std::chrono::system_clock::now();
auto t = std::chrono::system_clock::to_time_t(now);
std::tm tm{};
#ifdef _WIN32
gmtime_s(&tm, &t);
#else
gmtime_r(&t, &tm);
#endif
char buf[32];
std::strftime(buf, sizeof(buf), "%Y-%m-%d", &tm);
return buf;
}
/// Format a double with given decimal places; returns "?" for non-finite.
/// snprintf rather than an ostringstream: the reflection loop below calls this twelve times per
/// reflection, and building a stream (and its locale) per call dominated the time spent writing a
/// merged file. Same digits - both go through the C locale's %.*f.
std::string Fmt(double val, int decimals = 4) {
if (!std::isfinite(val))
return "?";
char buf[512];
const int n = std::snprintf(buf, sizeof(buf), "%.*f", decimals, val);
return std::string(buf, n);
}
/// Quote a CIF string value; returns "?" for empty.
std::string CifStr(const std::string& s) {
if (s.empty())
return "?";
// If it contains spaces or special chars, single-quote it
if (s.find(' ') != std::string::npos ||
s.find('\'') != std::string::npos ||
s.find('#') != std::string::npos)
return "'" + s + "'";
return s;
}
// One output row per reflection in the standard CCP4 anomalous layout: the merged mean (IMEAN / F)
// plus the two Bijvoet mates (I(+)/I(-), F(+)/F(-)). The merge keeps the two mates as separate rows
// (I+ under the ASU representative hkl, I- under -hkl); this collapses them into one row so the MTZ
// and SHELX writers share one row list. has_anom is set false when NO reflection carries an anomalous
// split (e.g. the stills path) - callers then omit the +/- columns.
struct MergedOutRow {
int h = 0, k = 0, l = 0;
float Imean = NAN, sImean = NAN, Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
float Fmean = NAN, sFmean = NAN, Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
int rfree = 0;
};
std::vector<MergedOutRow> BuildMergedRows(const std::vector<MergedReflection> &reflections,
const DiffractionExperiment &experiment,
bool &has_anom) {
std::vector<MergedOutRow> out_rows;
has_anom = true;
if (experiment.GetScalingSettings().GetMergeFriedel()) {
// Friedel-merged: IMEAN is the already-merged intensity (r.I). I(+)/I(-) are carried verbatim
// from the Bijvoet split the merge kept (rotation always does; scaled non-anomalously), so a weak
// anomalous signal is preserved without reprocessing. A reflection with only one mate, or a
// centric, gets a missing value (NaN) for the absent hand. When NO reflection has an anomalous
// split (e.g. the stills path, which does not compute one) the anomalous columns are omitted.
has_anom = std::any_of(reflections.begin(), reflections.end(),
[](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); });
out_rows.reserve(reflections.size());
for (const auto& r : reflections)
out_rows.push_back({r.h, r.k, r.l, r.I, r.sigma, r.I_plus, r.sigma_plus, r.I_minus,
r.sigma_minus, r.F, r.sigmaF, r.F_plus, r.sigmaF_plus, r.F_minus,
r.sigmaF_minus, r.rfree_flag ? 1 : 0});
} else {
// Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form
// IMEAN / F as their inverse-variance Friedel mean. A single generator gives both the group key
// (its hkl, identical for +hkl and -hkl) and which mate this row is (.plus).
const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupNumber().value_or(1));
struct AnomRow {
int h = 0, k = 0, l = 0;
float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
float Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
int rfree = 0;
};
std::map<std::tuple<int, int, int>, AnomRow> rows;
for (const auto& r : reflections) {
const HKLKey key = key_gen(r);
AnomRow& row = rows[{key.h, key.k, key.l}];
row.h = key.h; row.k = key.k; row.l = key.l;
row.rfree = r.rfree_flag ? 1 : 0;
if (key.plus) { row.Ip = r.I; row.sIp = r.sigma; row.Fp = r.F; row.sFp = r.sigmaF; }
else { row.Im = r.I; row.sIm = r.sigma; row.Fm = r.F; row.sFm = r.sigmaF; }
}
// Friedel-mean of the two mates by inverse variance (the single mate, if only one was measured).
const auto combine = [](float a, float sa, float b, float sb, float& val, float& sig) {
const bool ok_a = std::isfinite(a) && sa > 0.0f;
const bool ok_b = std::isfinite(b) && sb > 0.0f;
if (ok_a && ok_b) {
const double wa = 1.0 / (static_cast<double>(sa) * sa);
const double wb = 1.0 / (static_cast<double>(sb) * sb);
val = static_cast<float>((wa * a + wb * b) / (wa + wb));
sig = static_cast<float>(1.0 / std::sqrt(wa + wb));
} else if (ok_a) { val = a; sig = sa; }
else if (ok_b) { val = b; sig = sb; }
else { val = NAN; sig = NAN; }
};
out_rows.reserve(rows.size());
for (const auto& [hkl, row] : rows) {
float i_mean, sig_i_mean, f_mean, sig_f_mean;
combine(row.Ip, row.sIp, row.Im, row.sIm, i_mean, sig_i_mean);
combine(row.Fp, row.sFp, row.Fm, row.sFm, f_mean, sig_f_mean);
out_rows.push_back({row.h, row.k, row.l, i_mean, sig_i_mean, row.Ip, row.sIp, row.Im,
row.sIm, f_mean, sig_f_mean, row.Fp, row.sFp, row.Fm, row.sFm, row.rfree});
}
}
return out_rows;
}
} // namespace
void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
const UnitCell &unitCell,
const DiffractionExperiment &experiment,
const MergeStatistics &statistics,
const ErrorModelReport &error_model,
const TwinningAnalysisResult &twinning,
const std::string &filename,
size_t nthreads) {
std::ofstream out(filename);
if (!out)
throw std::runtime_error("WriteMmcifReflections: cannot open " + filename);
out << std::fixed;
// ---------- data block ----------
out << "data_sample" << "\n";
out << "#\n";
// ---------- _audit ----------
out << "_audit.revision_id 1\n";
out << "_audit.creation_date " << CurrentDateISO() << "\n";
out << "_audit.update_record 'Initial release'\n";
out << "#\n";
// ---------- _software ----------
out << "_software.name 'Rugnux'\n";
out << "_software.version " << CifStr(jfjoch_version()) << "\n";
out << "_software.classification 'data reduction'\n";
out << "#\n";
// ---------- _cell ----------
out << "_cell.length_a " << Fmt(unitCell.a, 3) << "\n";
out << "_cell.length_b " << Fmt(unitCell.b, 3) << "\n";
out << "_cell.length_c " << Fmt(unitCell.c, 3) << "\n";
out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 2) << "\n";
out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n";
out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\n";
auto *sg = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1));
if (sg == nullptr)
throw std::runtime_error("WriteMmcifReflections: invalid space group number");
// ---------- _symmetry ----------
out << "_symmetry.space_group_name_H-M " << CifStr(sg->hm) << "\n";
out << "_symmetry.Int_Tables_number " << sg->number << "\n";
out << "#\n";
// ---------- _diffrn_source / _diffrn_detector ----------
if (!experiment.GetSourceName().empty())
out << "_diffrn_source.pdbx_synchrotron_site " << CifStr(experiment.GetSourceName()) << "\n";
if (!experiment.GetInstrumentName().empty())
out << "_diffrn_source.pdbx_synchrotron_beamline " << CifStr(experiment.GetInstrumentName()) << "\n";
out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n";
out << "_diffrn_detector.detector " << CifStr(experiment.GetDetectorDescription()) << "\n";
// Detector geometry actually used for integration (refined, when geometry refinement ran - the rotation
// two-pass or the stills global refinement update it on experiment_ before the written pass). jfjoch_
// local-data-name items: feedback of the distance / beam centre the data was reduced with.
out << "_diffrn_detector.jfjoch_distance_mm " << Fmt(experiment.GetDetectorDistance_mm(), 4) << "\n";
out << "_diffrn_detector.jfjoch_beam_center_x_pxl " << Fmt(experiment.GetBeamX_pxl(), 2) << "\n";
out << "_diffrn_detector.jfjoch_beam_center_y_pxl " << Fmt(experiment.GetBeamY_pxl(), 2) << "\n";
out << "#\n";
// ---------- merging statistics (_reflns overall + _reflns_shell loop) ----------
// cc_half and r_meas are stored as fractions (0-1), which is the mmCIF convention. ISa (the
// Diederichs asymptotic I/sigma, 1/b of the a*sigma^2 + (b*I)^2 error model) and the twinning
// indicators below have no standard mmCIF item. They are written under the "jfjoch" reserved
// prefix (_reflns.jfjoch_*), the IUCr-sanctioned local-data-name extension for private items -
// NOT the "pdbx_" prefix, which is owned by the wwPDB PDBx/mmCIF dictionary and must not label
// items that dictionary does not define. (The other pdbx_ items here are genuine PDBx items.)
const auto mult = [](const MergeStatisticsShell &s) {
return s.unique_reflections > 0 ? static_cast<double>(s.total_observations) / s.unique_reflections : 0.0; };
const auto compl_pct = [](const MergeStatisticsShell &s) {
return s.possible_unique_reflections > 0
? 100.0 * static_cast<double>(s.unique_reflections) / s.possible_unique_reflections : 0.0; };
if (!statistics.shells.empty()) {
const auto &ov = statistics.overall;
// Anomalous signal-to-noise (SigAno) is written only when an anomalous split was made, so a
// non-anomalous merge keeps its previous stats block / shell-loop columns unchanged.
const bool has_anom = std::isfinite(ov.abs_diff_over_sigma_anomalous);
out << "_reflns.d_resolution_high " << Fmt(ov.d_min, 2) << "\n";
out << "_reflns.d_resolution_low " << Fmt(ov.d_max, 2) << "\n";
out << "_reflns.number_obs " << ov.unique_reflections << "\n";
out << "_reflns.pdbx_number_measured_all " << ov.total_observations << "\n";
out << "_reflns.pdbx_redundancy " << Fmt(mult(ov), 2) << "\n";
out << "_reflns.percent_possible_obs " << Fmt(compl_pct(ov), 1) << "\n";
out << "_reflns.pdbx_netI_over_sigmaI " << Fmt(ov.mean_i_over_sigma, 2) << "\n";
out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n";
out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n";
if (has_anom)
out << "_reflns.pdbx_absDiff_over_sigma_anomalous " << Fmt(ov.abs_diff_over_sigma_anomalous, 3)
<< " # SigAno = <|dano|>/<sigma(dano)>\n";
// ISa in XDS's sense: the whole-range 1/sqrt(a*b) of the error model below, so this item can
// be read straight against a CORRECT.LP. The strong-reflection asymptote - a tier XDS does not
// have, and always the more optimistic of the two - is written separately rather than here.
out << "_reflns.jfjoch_diffrn_ISa " << CifStr(error_model.isa)
<< " # 1/sqrt(a*b), the XDS convention\n";
if (!error_model.isa_asymptotic.empty())
out << "_reflns.jfjoch_diffrn_ISa_asymptotic " << CifStr(error_model.isa_asymptotic)
<< " # strong-reflection asymptote (Diederichs); rotation path only\n";
if (!error_model.a.empty())
out << "_reflns.jfjoch_error_model_a " << CifStr(error_model.a)
<< " # sigma^2 = a*(sigma0^2 + b*I^2), XDS convention\n";
if (!error_model.b.empty())
out << "_reflns.jfjoch_error_model_b " << CifStr(error_model.b) << "\n";
// Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's
// "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded.
if (std::isfinite(statistics.wilson_b) && statistics.wilson_b > 0.0)
out << "_reflns.B_iso_Wilson_estimate " << Fmt(statistics.wilson_b, 2) << "\n";
// Twinning indicators (no standard mmCIF item; same jfjoch local prefix as ISa above).
if (twinning.l_test_pairs > 0) {
out << "_reflns.jfjoch_L_test_mean_abs_L " << Fmt(twinning.mean_abs_l, 3)
<< " # Padilla-Yeates <|L|> (untwinned 0.500, perfect twin 0.375)\n";
out << "_reflns.jfjoch_L_test_mean_L_squared " << Fmt(twinning.mean_l_squared, 3)
<< " # <L^2> (untwinned 0.333, perfect twin 0.200)\n";
}
if (twinning.moment_reflections > 0)
out << "_reflns.jfjoch_second_moment_I " << Fmt(twinning.second_moment, 3)
<< " # <I^2>/<I>^2 (untwinned 2.00, perfect twin 1.50)\n";
// Radiation-damage monitor (rotation): the relative Debye-Waller B change from the first to the last
// frame (A^2). A large magnitude flags a dose-dependent resolution-scale change = radiation damage;
// positive is the typical direction (high-resolution intensity fades with dose). No standard mmCIF item.
if (std::isfinite(statistics.radiation_damage_delta_b))
out << "_reflns.jfjoch_radiation_damage_relative_B " << Fmt(statistics.radiation_damage_delta_b, 2)
<< " # relative-B first->last over the run (A^2); + = high-res fades with dose\n";
// Diffraction anisotropy. The eigen-decomposition of the anisotropy tensor has standard PDBx
// items; the eigenvalues there must be non-negative, and only the deviatoric part of the tensor
// is determined at all (its isotropic part is degenerate with the overall scale), so they are
// written relative to the weakest direction - eigenvalue_3 is 0 by construction and
// eigenvalue_1 is the anisotropic deltaB. The eigenvectors are in the PDB orthogonalisation
// convention, which is the one gemmi (and hence rugnux) uses throughout.
//
// The directional diffraction LIMITS are deliberately NOT written as
// _reflns.pdbx_aniso_diffraction_limit_*: the dictionary defines those as the semi-axes of an
// ellipsoid fitted to a diffraction cut-off surface, which is a different construction from the
// one below and one rugnux does not perform - it cuts nothing on a directional criterion. They
// go under the jfjoch local prefix with their own definition instead.
const auto &an = statistics.anisotropy;
if (an.n_reflections > 0 && an.n_cells > 0 && std::isfinite(an.delta_b)) {
out << "_reflns.pdbx_orthogonalization_convention pdb\n";
for (int i = 0; i < 3; ++i) {
out << "_reflns.pdbx_aniso_B_tensor_eigenvalue_" << (i + 1) << " "
<< Fmt(an.eigenvalue[i] - an.eigenvalue[2], 2) << "\n";
for (int j = 0; j < 3; ++j)
out << "_reflns.pdbx_aniso_B_tensor_eigenvector_" << (i + 1) << "_ortho[" << (j + 1)
<< "] " << Fmt(an.eigenvector[i][j], 4) << "\n";
}
out << "_reflns.jfjoch_aniso_delta_B " << Fmt(an.delta_b, 2)
<< " # range of the principal components (A^2), fitted on intensities\n";
if (std::isfinite(an.delta_b_linear))
out << "_reflns.jfjoch_aniso_delta_B_linear " << Fmt(an.delta_b_linear, 2)
<< " # deltaB implied by the s^2 slope alone: what a Debye-Waller B accounts for\n";
out << "_reflns.jfjoch_aniso_shape " << AnisotropyShapeCode(an.shape)
<< " # resolution signature of the directional deficit\n";
if (std::isfinite(an.floor))
out << "_reflns.jfjoch_aniso_floor " << Fmt(an.floor, 3)
<< " # deltaB this data set's own systematic error could manufacture (A^2)\n";
if (std::isfinite(an.significance))
out << "_reflns.jfjoch_aniso_significance " << Fmt(an.significance, 2)
<< " # deltaB(linear) / floor\n";
out << "_reflns.jfjoch_aniso_verdict " << AnisotropyVerdictCode(an.verdict)
<< " # DETECTED / NOT_DETECTED / CANNOT_DETERMINE, at deltaB(linear)/floor > 3.5\n";
for (int i = 0; i < 3; ++i)
if (std::isfinite(an.d_min_axis[i]))
out << "_reflns.jfjoch_aniso_d_min_" << (i + 1) << " "
<< Fmt(an.d_min_axis[i], 2)
<< " # <I/sigma(I)> = 2 in a 20 deg cone about eigenvector " << (i + 1)
<< (an.d_min_censored[i] ? "; at the edge of the measured data, so a bound" : "")
<< "\n";
}
out << "#\n";
// Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B
// per rotation-range batch, measured before any correction. rotation_start_deg = id * batch_deg.
if (!statistics.radiation_damage_b_batch.empty()) {
out << "loop_\n";
out << "_jfjoch_radiation_damage_batch.id\n";
out << "_jfjoch_radiation_damage_batch.rotation_start_deg\n";
out << "_jfjoch_radiation_damage_batch.relative_B\n";
for (size_t i = 0; i < statistics.radiation_damage_b_batch.size(); ++i)
out << " " << (i + 1) << " "
<< Fmt(static_cast<double>(i) * statistics.radiation_damage_batch_deg, 1) << " "
<< Fmt(statistics.radiation_damage_b_batch[i], 2) << "\n";
out << "#\n";
}
out << "loop_\n";
out << "_reflns_shell.d_res_high\n";
out << "_reflns_shell.d_res_low\n";
out << "_reflns_shell.number_measured_obs\n";
out << "_reflns_shell.number_unique_obs\n";
out << "_reflns_shell.pdbx_redundancy\n";
out << "_reflns_shell.percent_possible_obs\n";
out << "_reflns_shell.meanI_over_sigI_obs\n";
out << "_reflns_shell.pdbx_Rrim_I_all\n";
out << "_reflns_shell.pdbx_CC_half\n";
if (has_anom)
out << "_reflns_shell.pdbx_absDiff_over_sigma_anomalous\n";
for (const auto &s : statistics.shells) {
if (s.unique_reflections == 0)
continue;
out << Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " "
<< s.total_observations << " " << s.unique_reflections << " "
<< Fmt(mult(s), 2) << " " << Fmt(compl_pct(s), 1) << " "
<< Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_meas, 4) << " " << Fmt(s.cc_half, 4);
if (has_anom)
out << " " << Fmt(s.abs_diff_over_sigma_anomalous, 3);
out << "\n";
}
out << "#\n";
}
// ---------- _refln loop ----------
out << "loop_\n";
out << "_refln.index_h\n";
out << "_refln.index_k\n";
out << "_refln.index_l\n";
out << "_refln.intensity_meas\n";
out << "_refln.intensity_sigma\n";
out << "_refln.pdbx_I_plus\n";
out << "_refln.pdbx_I_plus_sigma\n";
out << "_refln.pdbx_I_minus\n";
out << "_refln.pdbx_I_minus_sigma\n";
out << "_refln.F_meas_au\n";
out << "_refln.F_meas_sigma_au\n";
out << "_refln.pdbx_F_plus\n";
out << "_refln.pdbx_F_plus_sigma\n";
out << "_refln.pdbx_F_minus\n";
out << "_refln.pdbx_F_minus_sigma\n";
out << "_refln.status_free\n";
out << "_refln.status\n";
// One row per unique reflection, twelve formatted floats each - tens of megabytes on a crowded
// crystal, and the largest single-threaded stretch left in a run. Nothing about a row depends on
// any other, so each worker formats its own block into its own string and the blocks go to the
// file in order. The columns are written exactly as the stream wrote them: the same "%.4f" (or
// "?" where the value is not finite), right-aligned in the same width.
{
const size_t nrow = reflections.size();
const size_t nw = std::max<size_t>(nthreads, 1);
const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
std::vector<std::string> block(nch);
ParallelChunks(nch, nw, [&](int tlo, int thi) {
for (int t = tlo; t < thi; ++t) {
const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch;
std::string &s = block[t];
s.reserve((hi - lo) * 208);
const auto column = [&s](const std::string &v, size_t width) {
if (v.size() < width) s.append(width - v.size(), ' ');
s.append(v);
s.push_back(' ');
};
for (size_t i = lo; i < hi; ++i) {
const auto &r = reflections[i];
column(std::to_string(r.h), 5);
column(std::to_string(r.k), 5);
column(std::to_string(r.l), 5);
column(Fmt(r.I, 4), 14);
column(Fmt(r.sigma, 4), 14);
column(Fmt(r.I_plus, 4), 14);
column(Fmt(r.sigma_plus, 4), 14);
column(Fmt(r.I_minus, 4), 14);
column(Fmt(r.sigma_minus, 4), 14);
column(Fmt(r.F, 4), 14);
column(Fmt(r.sigmaF, 4), 14);
column(Fmt(r.F_plus, 4), 14);
column(Fmt(r.sigmaF_plus, 4), 14);
column(Fmt(r.F_minus, 4), 14);
column(Fmt(r.sigmaF_minus, 4), 14);
s.push_back(r.rfree_flag ? '1' : '0');
s.append(" o\n"); // 'o' = observed
}
}
});
for (const std::string &s : block)
out.write(s.data(), static_cast<std::streamsize>(s.size()));
}
out << "#\n";
out << "# End of reflections\n";
out.close();
}
void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
const UnitCell &unitCell,
const DiffractionExperiment &experiment,
const std::string &filename) {
gemmi::Mtz mtz;
// Optional but recommended metadata
mtz.spacegroup = gemmi::find_spacegroup_by_number(
experiment.GetSpaceGroupNumber().value_or(1));
mtz.set_cell_for_all(unitCell);
// Producing-software provenance in the MTZ header (title + HISTORY, the CCP4 convention).
mtz.title = "Rugnux merged reflections";
mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction");
// Add dataset
gemmi::Mtz::Dataset& ds = mtz.add_dataset("native");
ds.crystal_name = experiment.GetSampleName();
ds.wavelength = experiment.GetWavelength_A();
const int dataset_id = ds.id;
// One row per reflection in the CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for
// the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly.
bool has_anom = true;
const std::vector<MergedOutRow> out_rows = BuildMergedRows(reflections, experiment, has_anom);
mtz.add_column("H", 'H', dataset_id, -1, false);
mtz.add_column("K", 'H', dataset_id, -1, false);
mtz.add_column("L", 'H', dataset_id, -1, false);
mtz.add_column("IMEAN", 'J', dataset_id, -1, false);
mtz.add_column("SIGIMEAN", 'Q', dataset_id, -1, false);
if (has_anom) {
mtz.add_column("I(+)", 'K', dataset_id, -1, false);
mtz.add_column("SIGI(+)", 'M', dataset_id, -1, false);
mtz.add_column("I(-)", 'K', dataset_id, -1, false);
mtz.add_column("SIGI(-)", 'M', dataset_id, -1, false);
}
mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude
mtz.add_column("SIGF", 'Q', dataset_id, -1, false);
if (has_anom) {
mtz.add_column("F(+)", 'G', dataset_id, -1, false);
mtz.add_column("SIGF(+)", 'L', dataset_id, -1, false);
mtz.add_column("F(-)", 'G', dataset_id, -1, false);
mtz.add_column("SIGF(-)", 'L', dataset_id, -1, false);
}
mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false);
mtz.nreflections = static_cast<int>(out_rows.size());
mtz.data.reserve(out_rows.size() * (has_anom ? 16 : 8));
for (const auto& row : out_rows) {
mtz.data.push_back(static_cast<float>(row.h));
mtz.data.push_back(static_cast<float>(row.k));
mtz.data.push_back(static_cast<float>(row.l));
mtz.data.push_back(row.Imean);
mtz.data.push_back(row.sImean);
if (has_anom) {
mtz.data.push_back(row.Ip);
mtz.data.push_back(row.sIp);
mtz.data.push_back(row.Im);
mtz.data.push_back(row.sIm);
}
mtz.data.push_back(row.Fmean);
mtz.data.push_back(row.sFmean);
if (has_anom) {
mtz.data.push_back(row.Fp);
mtz.data.push_back(row.sFp);
mtz.data.push_back(row.Fm);
mtz.data.push_back(row.sFm);
}
mtz.data.push_back(static_cast<float>(row.rfree));
}
mtz.write_to_file(filename);
}
void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
const DiffractionExperiment &experiment,
const std::string &filename,
size_t nthreads) {
bool has_anom = true;
const std::vector<MergedOutRow> rows = BuildMergedRows(reflections, experiment, has_anom);
// SHELX HKLF 4 (SHELXC / ANODE input): fixed FORMAT(3I4,2F8.2), one record per reflection as
// h k l I sigma(I). The Bijvoet mates are written separately - I(+) at +hkl, I(-) at -hkl - so the
// anomalous differences survive; a reflection with no anomalous split is written once as its mean.
// Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale
// is irrelevant to SHELXC / ANODE, which use only ratios); I and sigma share the scale, so the
// anomalous signal is untouched. The file ends with a 0 0 0 terminator record.
const auto usable = [](float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; };
double max_abs = 0.0;
for (const auto& r : rows) {
if (usable(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)});
if (usable(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)});
if (!usable(r.Ip, r.sIp) && !usable(r.Im, r.sIm) && usable(r.Imean, r.sImean))
max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)});
}
const double scale = (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0;
std::ofstream out(filename);
if (!out)
throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename);
// Up to two records per reflection, five formatted numbers each. Built in parallel into per-worker
// blocks and handed to the file in order, exactly as the mmCIF rows are; "%.2f" right-aligned in
// the fixed field is what `fixed` + `setprecision(2)` + `setw` made the stream write.
const auto column = [](std::string &s, const std::string &v, size_t width) {
if (v.size() < width) s.append(width - v.size(), ' ');
s.append(v);
};
const auto num2 = [](double v) {
char b[64];
const int n = std::snprintf(b, sizeof b, "%.2f", v);
return std::string(b, static_cast<size_t>(std::clamp(n, 0, static_cast<int>(sizeof b) - 1)));
};
const auto emit = [&column, &num2, scale](std::string &s, int h, int k, int l, float I, float sigma) {
column(s, std::to_string(h), 4);
column(s, std::to_string(k), 4);
column(s, std::to_string(l), 4);
column(s, num2(scale * I), 8);
column(s, num2(scale * sigma), 8);
s.push_back('\n');
};
{
const size_t nrow = rows.size();
const size_t nw = std::max<size_t>(nthreads, 1);
const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
std::vector<std::string> block(nch);
ParallelChunks(nch, nw, [&](int tlo, int thi) {
for (int t = tlo; t < thi; ++t) {
const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch;
std::string &s = block[t];
s.reserve((hi - lo) * 2 * 29);
for (size_t i = lo; i < hi; ++i) {
const auto &r = rows[i];
const bool plus = usable(r.Ip, r.sIp);
const bool minus = usable(r.Im, r.sIm);
if (plus) emit(s, r.h, r.k, r.l, r.Ip, r.sIp);
if (minus) emit(s, -r.h, -r.k, -r.l, r.Im, r.sIm);
if (!plus && !minus && usable(r.Imean, r.sImean))
emit(s, r.h, r.k, r.l, r.Imean, r.sImean);
}
}
});
for (const std::string &s : block)
out.write(s.data(), static_cast<std::streamsize>(s.size()));
}
std::string tail;
emit(tail, 0, 0, 0, 0.0f, 0.0f); // HKLF-4 end-of-data marker
out.write(tail.data(), static_cast<std::streamsize>(tail.size()));
out.close();
}
namespace {
// The detector position an observation is written at: a reflection too weak to have a measured
// centroid still has a predicted one.
float DetectorX(const Reflection &r) { return std::isfinite(r.observed_x) ? r.observed_x : r.predicted_x; }
float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.observed_y : r.predicted_y; }
// Sum each rocking event into one full observation. A rotation reflection is integrated image by
// image, so it arrives here as a run of partials over consecutive frames; the run is cut where the
// 3D combine cuts it - same raw hkl, frames no further apart than MAX_FRAME_GAP - so the exported
// file and rugnux's own merge see exactly the same events.
// The parts are added, plainly, with their variances in quadrature, which is what every other
// rotation program writes as a full. Nothing is divided by the partiality: FRACTIONCALC carries the
// summed rocking-curve fraction, and an event the sweep cut short says so there rather than being
// scaled up or dropped here.
// The metadata a full carries is the partiality-weighted mean of its parts - the centroid of the
// rocking curve, which is what XD/YD/ZD mean in an XDS file and what a scale model downstream reads
// these columns for. An intensity weighting is the other candidate and it does not survive weak
// data: the parts of a weak reflection scatter about zero, so on the crystals this was measured on
// 8-41% of events have parts summing to zero or less and 10-21% of the intensity-weighted centroids
// fall outside the event's own frame range. Partialities are the rocking curve's own weights and are
// positive by construction.
// min_partiality is the combine's own floor on the assembled full (--min-partiality, default 0.02,
// 0 = off). An event that caught a thousandth of its rocking curve is not a measurement of that
// reflection - rugnux does not merge one either - and writing it as a full hands the reading program
// a whole observation whose intensity is noise and whose sigma is small, so it is weighted heavily.
// On a crystal whose rocking curves span twenty frames those events are 7% of the file and cost
// POINTLESS the point group. Everything above the floor is written with its honest FRACTIONCALC,
// truncated edge-of-sweep events included, for the reader to cut where it wants.
std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &outcomes,
double min_partiality) {
constexpr float MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge's: what makes one rocking event
// The sort key travels with the part instead of being read back through the pointer, the way the
// merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts
// and an indirect compare is a cache miss on every one of them. The keys are the same values in
// the same order, so introsort makes the same comparisons and the same swaps and leaves the same
// order - which matters, because two parts can genuinely share (h,k,l) and image_number and the
// event sums below are floating point.
struct Part {
int32_t h, k, l;
float image_number;
const Reflection *r;
};
std::vector<Part> parts;
size_t n_parts = 0;
for (const auto &outcome : outcomes)
n_parts += outcome.reflections.size();
parts.reserve(n_parts);
for (const auto &outcome : outcomes)
for (const auto &r : outcome.reflections)
parts.push_back({r.h, r.k, r.l, r.image_number, &r});
std::sort(parts.begin(), parts.end(), [](const Part &a, const Part &b) {
return std::tie(a.h, a.k, a.l, a.image_number) < std::tie(b.h, b.k, b.l, b.image_number);
});
std::vector<Reflection> fulls;
for (size_t i = 0; i < parts.size(); ) {
size_t j = i + 1;
while (j < parts.size() && parts[j].h == parts[i].h && parts[j].k == parts[i].k
&& parts[j].l == parts[i].l
&& parts[j].image_number - parts[j - 1].image_number <= MAX_FRAME_GAP)
++j;
double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, sum_var_bkg = 0.0;
double p_rlp = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0, p_delta_phi = 0.0, p_zeta = 0.0,
p_bkg = 0.0;
for (size_t m = i; m < j; ++m) {
const Reflection &r = *parts[m].r;
const double p = r.partiality;
sum_p += p;
sum_I += static_cast<double>(r.I) * r.rlp;
sum_var += static_cast<double>(r.sigma) * r.sigma * r.rlp * r.rlp;
sum_var_bkg += static_cast<double>(r.var_bkg) * r.rlp * r.rlp;
p_rlp += p * r.rlp;
p_frame += p * r.image_number;
p_x += p * DetectorX(r);
p_y += p * DetectorY(r);
p_delta_phi += p * r.delta_phi_deg;
p_zeta += p * r.zeta;
p_bkg += p * r.bkg;
}
Reflection full = *parts[i].r;
i = j;
if (sum_p < min_partiality)
continue;
// The Lorentz-polarization factor is applied by the writer, which multiplies I by rlp, so
// divide the event's own factor back out of the sums here. LP is the same geometry for every
// part of one event to a median 2e-4, so the file's I/LP is still the raw count sum.
full.rlp = static_cast<float>(p_rlp / sum_p);
full.I = static_cast<float>(sum_I / full.rlp);
full.sigma = static_cast<float>(std::sqrt(sum_var) / full.rlp);
full.var_bkg = static_cast<float>(sum_var_bkg / (static_cast<double>(full.rlp) * full.rlp));
full.partiality = static_cast<float>(sum_p);
full.image_number = static_cast<float>(p_frame / sum_p);
full.observed_x = static_cast<float>(p_x / sum_p);
full.observed_y = static_cast<float>(p_y / sum_p);
full.delta_phi_deg = static_cast<float>(p_delta_phi / sum_p);
full.zeta = static_cast<float>(p_zeta / sum_p);
full.bkg = static_cast<float>(p_bkg / sum_p);
fulls.push_back(full);
}
return fulls;
}
} // namespace
void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes,
const UnitCell &unitCell,
const DiffractionExperiment &experiment,
bool sum_partials,
const std::string &filename) {
gemmi::Mtz mtz;
mtz.spacegroup = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1));
mtz.set_cell_for_all(unitCell);
mtz.title = "Rugnux unmerged observations";
mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction");
mtz.add_base(); // the HKL_base dataset and the H K L columns
gemmi::Mtz::Dataset &ds = mtz.add_dataset("unmerged");
ds.crystal_name = experiment.GetSampleName();
ds.wavelength = experiment.GetWavelength_A();
// Every data column belongs to HKL_base and only the batches carry a dataset id, as in a
// POINTLESS file. Down to FLAG this is POINTLESS's own column set; the four after it are what
// rugnux measures beyond it - the offset of the reflection from the centre of its rocking curve,
// the Lorentz geometry of that curve, and the background that was subtracted.
mtz.add_column("M/ISYM", 'Y', 0, -1, false);
mtz.add_column("BATCH", 'B', 0, -1, false);
mtz.add_column("I", 'J', 0, -1, false);
mtz.add_column("SIGI", 'Q', 0, -1, false);
mtz.add_column("FRACTIONCALC", 'R', 0, -1, false);
mtz.add_column("XDET", 'R', 0, -1, false);
mtz.add_column("YDET", 'R', 0, -1, false);
mtz.add_column("ROT", 'R', 0, -1, false);
mtz.add_column("LP", 'R', 0, -1, false);
mtz.add_column("FLAG", 'I', 0, -1, false);
mtz.add_column("DELPHI", 'R', 0, -1, false);
mtz.add_column("ZETA", 'R', 0, -1, false);
mtz.add_column("BGMEAN", 'R', 0, -1, false);
mtz.add_column("BGVAR", 'R', 0, -1, false);
const auto gon = experiment.GetGoniometer();
// A scanning goniometer means rugnux integrated a rocking curve image by image, so each
// observation is one PART of a reflection unless the parts are summed here; on a still each is a
// whole measurement either way, and consecutive stills are different crystals, so there is
// nothing to sum there. Which of the two the file holds has to be declared in two places, because
// that is where POINTLESS and AIMLESS look for it: the batch header's data type, and the M flag
// that M/ISYM packs as 256*M + ISYM. Told nothing, they take each partial for a whole reflection
// and neither the symmetry determination nor the scaling survives it.
const bool scanning = gon && gon->IsScanning();
const bool partials = scanning && !sum_partials;
const float wedge_deg = gon ? gon->GetWedge_deg() : 0.0f;
const auto phi_start_deg = [&](float image_number) {
return gon ? gon->GetAngle_deg(image_number) : 0.0f;
};
// H K L are the ASU indices and M/ISYM says which symmetry operation (and which Friedel hand) got
// them there, so the index the reflection was actually measured at is recoverable - that is the
// crystal-frame information careless scales on, and what makes the file unmerged rather than a
// merge waiting to happen.
// I and SIGI are the integrated intensity with the Lorentz-polarization factor applied and
// nothing else, which is what IOBS means in every unmerged format (LP records the factor, so the
// raw counts are I/LP). LP is geometry, not a scale, and a program that reads this file has no
// way to recover it. The partiality is NOT divided out - that is a scale, FRACTIONCALC carries
// it, and every program this file is for wants to handle it its own way.
gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
std::set<int> batch_numbers;
// Lattice-centring absences are integrated on purpose - prediction runs in P so the space-group
// search can confirm or disprove the centring - but they are not written here. POINTLESS reads
// the interstitial nodes, takes the lattice for primitive and demotes the space group, and they
// roughly halve AIMLESS's <I/sigma>. Screw and glide absences ARE written: they are the evidence
// the space group was chosen on, so a program reading this file can check that choice, and
// removing them turns its test into an assumption. XDS and DIALS draw the line in the same place.
const char centering = mtz.spacegroup ? mtz.spacegroup->hm[0] : 'P';
const auto add_row = [&](const Reflection &r) {
if (systematic_absence(r.h, r.k, r.l, centering))
return;
std::array<int, 3> hkl{r.h, r.k, r.l};
const int isym = hkl_mover.move_to_asu(hkl);
// A summed full's image_number is its rocking-curve centroid, so this is the batch the
// reflection is centred on - which is what a batch means for a full everywhere else.
const int batch = 1 + static_cast<int>(std::lround(r.image_number));
batch_numbers.insert(batch);
mtz.data.push_back(static_cast<float>(hkl[0]));
mtz.data.push_back(static_cast<float>(hkl[1]));
mtz.data.push_back(static_cast<float>(hkl[2]));
mtz.data.push_back(static_cast<float>((partials ? 256 : 0) + isym));
mtz.data.push_back(static_cast<float>(batch));
mtz.data.push_back(r.I * r.rlp);
mtz.data.push_back(r.sigma * r.rlp);
mtz.data.push_back(r.partiality);
mtz.data.push_back(DetectorX(r));
mtz.data.push_back(DetectorY(r));
mtz.data.push_back(phi_start_deg(r.image_number) + wedge_deg / 2.0f);
mtz.data.push_back(r.rlp);
mtz.data.push_back(0.0f); // FLAG: nothing here is a rejected observation
mtz.data.push_back(r.delta_phi_deg);
mtz.data.push_back(r.zeta);
mtz.data.push_back(r.bkg);
mtz.data.push_back(r.var_bkg);
};
if (scanning && sum_partials) {
for (const auto &r : SumRockingEvents(outcomes,
experiment.GetScalingSettings().GetMinPartiality()))
add_row(r);
} else {
for (const auto &outcome : outcomes)
for (const auto &r : outcome.reflections)
add_row(r);
}
mtz.nreflections = static_cast<int>(mtz.data.size() / mtz.columns.size());
// The batch header's orientation matrix is the crystal at rotation angle zero - each batch's own
// PHISTT is applied on top of it - but the lattice stored with an outcome is the crystal as it
// stood on that image. Turn the first indexed one back by its own angle to get the orientation of
// the sweep, which is the one matrix POINTLESS also writes into every batch.
std::optional<CrystalLattice> lattice_at_zero;
std::optional<float> mosaicity_deg;
for (const auto &outcome : outcomes) {
if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f)
continue;
const float mid_deg = phi_start_deg(outcome.reflections.front().image_number) + wedge_deg / 2.0f;
lattice_at_zero = gon ? outcome.latt.Multiply(gon->GetTransformationAngle(mid_deg)) : outcome.latt;
mosaicity_deg = outcome.mosaicity_deg;
break;
}
// The batch header is written in the "Cambridge" frame - z along the principal rotation axis, x
// along the beam - while the jfjoch lab frame has the beam along +z, so the two are related by a
// rotation. These three lab-frame vectors are the Cambridge axes; a lab vector's components in
// that frame are its dot products with them. A still has no rotation axis, and any axis across
// the beam then defines the frame just as consistently.
// The axis is NEGATED: rugnux turns an observation made at angle phi back to phi = 0 by rotating
// it by +phi about the goniometer axis, so the crystal itself turns by -phi about it, while the
// MTZ batch header's scan axis is the one a batch's own increasing PHI turns the crystal about.
const Coord beam = experiment.GetDiffractionGeometry().GetScatteringVector().Normalize();
const Coord z_cam = gon ? -gon->GetAxis().Normalize() : Coord(0, 1, 0);
const Coord y_cam = (z_cam % beam).Normalize();
const Coord x_cam = (y_cam % z_cam).Normalize();
gemmi::Mtz::Batch batch;
batch.title = "Rugnux";
batch.axes.emplace_back("PHI");
batch.set_dataset_id(ds.id);
batch.ints[12] = 1; // ncryst
batch.ints[14] = partials ? 1 : 2; // ldtype: oscillation data (2D spots) / area detector (3D)
batch.ints[15] = 1; // jsaxs: the goniostat scan axis
batch.ints[17] = 1; // ngonax
batch.ints[19] = 1; // ndet
batch.set_cell(unitCell);
if (lattice_at_zero) {
// Orientation matrix U, built from the reciprocal axes and stored column by column in
// Cambridge components, as gemmi's XDS_ASCII converter builds it (gemmi/xds2mtz.hpp).
const Coord ar = lattice_at_zero->Astar().Normalize();
const Coord cr = (ar % lattice_at_zero->Bstar()).Normalize();
const Coord u[3] = {ar, cr % ar, cr};
for (int i = 0; i < 3; i++) {
batch.floats[6 + 3 * i] = u[i] * x_cam;
batch.floats[7 + 3 * i] = u[i] * y_cam;
batch.floats[8 + 3 * i] = u[i] * z_cam;
}
}
batch.floats[21] = mosaicity_deg.value_or(0.0f); // crydat(0), the reflecting range
batch.floats[40] = 1.0f; // scanax = [0, 0, 1]: the rotation axis IS z in the Cambridge frame
batch.floats[47] = wedge_deg;
batch.floats[61] = 1.0f; // e1 = scanax, the only goniostat axis
batch.floats[80] = -1.0f; // idealised source vector, antiparallel to the beam
batch.floats[83] = -(beam * x_cam); // s0, the source vector of the geometry as it really stands
batch.floats[84] = -(beam * y_cam);
batch.floats[85] = -(beam * z_cam);
batch.set_wavelength(experiment.GetWavelength_A());
batch.floats[111] = experiment.GetDetectorDistance_mm();
batch.floats[113] = 1.0f; // detector limits, in pixels
batch.floats[114] = static_cast<float>(experiment.GetXPixelsNum());
batch.floats[115] = 1.0f;
batch.floats[116] = static_cast<float>(experiment.GetYPixelsNum());
for (const int number : batch_numbers) {
batch.number = number;
batch.floats[36] = phi_start_deg(static_cast<float>(number - 1)); // phistt
batch.floats[37] = batch.floats[36] + wedge_deg; // phiend
mtz.batches.push_back(batch);
}
mtz.sort(5); // by H K L M/ISYM BATCH, the order POINTLESS leaves an unmerged file in
mtz.write_to_file(filename);
}
void WriteReflections(const std::vector<MergedReflection> &reflections,
const UnitCell &unitCell,
const DiffractionExperiment &experiment,
const MergeStatistics &statistics,
const ErrorModelReport &error_model,
const TwinningAnalysisResult &twinning,
const std::string &filename,
size_t nthreads) {
// Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 /
// phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 as
// the SHELXC / ANODE substructure-solution input).
WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz");
WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning,
filename + ".cif", nthreads);
WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads);
}
Reference in New Issue View Git Blame Copy Permalink