From bb5aaca68ea80853fb3363fea67a7dea1192c857 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 25 Aug 2026 16:18:14 +0200 Subject: [PATCH] rugnux: sum each rotation reflection's partials into one full on export --export-unmerged now writes one row per reflection with its partials summed, which is how a rotation dataset is normally handed to a scaling program. The previous behaviour is still available as --export-unmerged-partials. Partials made the file's usefulness depend on the reading program's partial-assembly settings. On an orthorhombic case only 31% of partial sets passed aimless's default PARTIALS test, so it merged 2% of the file: multiplicity 4.7 against our own 12.4, CC1/2 0.788, and CC against our merged intensities 0.806. Summed, aimless uses 89% of the file and rejects nothing: multiplicity 12.5, CC1/2 0.931, completeness 96.1 to 98.0, pointless Laue confidence 0.669 to 0.792, CC 0.848. The other two test crystals are ties overall and win their outer shells, and the anomalous signal comes out stronger from the fulls. iotbx.merging_statistics, which never assembled partials at all, now reports real merging statistics rather than statistics of parts. Events are grouped by Combine()'s own rule so the file and the merge assemble exactly the same reflections; the sum is plain with variances in quadrature and no partiality division, which reproduces the merge's full times FRACTIONCALC to a few percent. Per-part metadata is combined on a partiality-weighted centroid rather than an intensity-weighted one: 8-41% of events have parts summing to zero or less, where intensity weights are undefined or sign-flipped, and 10-21% of intensity-weighted centroids fall outside the event's own frame range. Events below --min-partiality are not written. Declaring a reflection whose rocking curve caught a thousandth of itself as a complete measurement gives a scaling program noise with a small sigma and it weights it heavily; there is no partial test to catch it once it is a full. This is the same constant already applied to the same quantity in the merge, and it is what takes the orthorhombic case from CC 0.219 to 0.848. --min-partiality 0 writes every event. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- docs/CHANGELOG.md | 2 +- image_analysis/WriteReflections.cpp | 170 ++++++++++++++++++++++------ image_analysis/WriteReflections.h | 21 ++-- rugnux/Rugnux.cpp | 19 +++- rugnux/Rugnux.h | 7 +- rugnux/rugnux_cli.cpp | 18 ++- 6 files changed, 186 insertions(+), 51 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b6fe74f3..e154e1dd 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,7 +3,7 @@ ### 1.0.0-rc.164 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. -* rugnux: `--export-unmerged ` writes the integrated observations as an unmerged MTZ - readable by aimless, pointless, careless and `iotbx.merging_statistics` - in `--mode mx` and `--mode scale` alike. +* rugnux: `--export-unmerged ` writes the integrated observations as an unmerged MTZ - readable by aimless, pointless, careless and `iotbx.merging_statistics` - in `--mode mx` and `--mode scale` alike, with each rotation reflection's partials summed into one full; `--export-unmerged-partials ` writes one row per image instead. ### 1.0.0-rc.163 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index cd72a274..5a402a21 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -538,9 +538,104 @@ void WriteShelxHklReflections(const std::vector &reflections, 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 SumRockingEvents(const std::vector &outcomes, + double min_partiality) { + constexpr float MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge's: what makes one rocking event + + std::vector parts; + for (const auto &outcome : outcomes) + for (const auto &r : outcome.reflections) + parts.push_back(&r); + std::sort(parts.begin(), parts.end(), [](const Reflection *a, const Reflection *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 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]; + const double p = r.partiality; + sum_p += p; + sum_I += static_cast(r.I) * r.rlp; + sum_var += static_cast(r.sigma) * r.sigma * r.rlp * r.rlp; + sum_var_bkg += static_cast(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]; + 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(p_rlp / sum_p); + full.I = static_cast(sum_I / full.rlp); + full.sigma = static_cast(std::sqrt(sum_var) / full.rlp); + full.var_bkg = static_cast(sum_var_bkg / (static_cast(full.rlp) * full.rlp)); + full.partiality = static_cast(sum_p); + full.image_number = static_cast(p_frame / sum_p); + full.observed_x = static_cast(p_x / sum_p); + full.observed_y = static_cast(p_y / sum_p); + full.delta_phi_deg = static_cast(p_delta_phi / sum_p); + full.zeta = static_cast(p_zeta / sum_p); + full.bkg = static_cast(p_bkg / sum_p); + fulls.push_back(full); + } + return fulls; +} + +} // namespace + void WriteUnmergedMtzReflections(const std::vector &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)); @@ -572,13 +667,15 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes 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 every - // observation here is one PART of a reflection; on a still each is a whole measurement. The - // difference 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 - // neither, they take each partial for a whole reflection and neither the symmetry determination - // nor the scaling survives it. - const bool partials = gon && gon->IsScanning(); + // 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; @@ -595,32 +692,39 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes // it, and every program this file is for wants to handle it its own way. gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup); std::set batch_numbers; - for (const auto &outcome : outcomes) { - for (const auto &r : outcome.reflections) { - std::array hkl{r.h, r.k, r.l}; - const int isym = hkl_mover.move_to_asu(hkl); - const int batch = 1 + static_cast(std::lround(r.image_number)); - batch_numbers.insert(batch); - mtz.data.push_back(static_cast(hkl[0])); - mtz.data.push_back(static_cast(hkl[1])); - mtz.data.push_back(static_cast(hkl[2])); - mtz.data.push_back(static_cast((partials ? 256 : 0) + isym)); - mtz.data.push_back(static_cast(batch)); - mtz.data.push_back(r.I * r.rlp); - mtz.data.push_back(r.sigma * r.rlp); - mtz.data.push_back(r.partiality); - // A reflection too weak to have a measured centroid still has a predicted position, and a - // detector coordinate is what the scale models downstream read this column for. - mtz.data.push_back(std::isfinite(r.observed_x) ? r.observed_x : r.predicted_x); - mtz.data.push_back(std::isfinite(r.observed_y) ? r.observed_y : r.predicted_y); - 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); - } + const auto add_row = [&](const Reflection &r) { + std::array 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(std::lround(r.image_number)); + batch_numbers.insert(batch); + mtz.data.push_back(static_cast(hkl[0])); + mtz.data.push_back(static_cast(hkl[1])); + mtz.data.push_back(static_cast(hkl[2])); + mtz.data.push_back(static_cast((partials ? 256 : 0) + isym)); + mtz.data.push_back(static_cast(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(mtz.data.size() / mtz.columns.size()); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 997b4d97..349813f3 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -47,17 +47,22 @@ void WriteShelxHklReflections(const std::vector &reflections, const std::string &filename, size_t nthreads); -// Unmerged observations, one row per integrated reflection, in the column and batch-header layout -// POINTLESS writes: aimless, pointless, careless and iotbx.merging_statistics all read that layout. -// H K L are the ASU indices and M/ISYM recovers the index the reflection was measured at (which is -// what careless needs to see the crystal frame) and says whether the observation is a partial. -// One observation per integrated box, so on a rotation run a reflection arrives as a run of partials -// over consecutive batches, for the reader to sum. The intensities carry the Lorentz-polarization -// factor and nothing else - the partiality and the per-image scale are left for the reading program -// to fit, since every program this file is for fits a scale model of its own. +// Unmerged observations in the column and batch-header layout POINTLESS writes: aimless, pointless, +// careless and iotbx.merging_statistics all read that layout. H K L are the ASU indices and M/ISYM +// recovers the index the reflection was measured at (which is what careless needs to see the crystal +// frame) and says whether the observation is a partial. +// sum_partials: add the partials of each rocking event into one full, written at the batch of the +// event's centroid with the summed rocking-curve fraction in FRACTIONCALC - the plain sum every +// rotation program writes, over the events the combine would assemble (--min-partiality). False +// writes one row per integrated box instead, flagged as partials for the reading program to sum. +// Stills have no rocking events and are unaffected either way. +// The intensities carry the Lorentz-polarization factor and nothing else - the partiality and the +// per-image scale are left for the reading program to fit, since every program this file is for fits +// a scale model of its own. void WriteUnmergedMtzReflections(const std::vector &outcomes, const UnitCell &unitCell, const DiffractionExperiment &experiment, + bool sum_partials, const std::string &filename); void WriteReflections(const std::vector &reflections, diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 00d827c1..815f0332 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2990,12 +2990,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // group is settled and the file has to be written in it, but from the integrated observations // rather than the merged ones: the partiality and the per-image scale are left for the reading // program, which fits a scale model of its own. - if (full && !cancelled_ && write_output && !geometry_prepass - && !config_.unmerged_export_path.empty() && result.consensus_cell) { - if (observer) observer->OnPhase("Writing unmerged reflections"); - WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, - experiment_, config_.unmerged_export_path); - logger.Info("Unmerged observations written to {}", config_.unmerged_export_path); + if (full && !cancelled_ && write_output && !geometry_prepass && result.consensus_cell) { + if (!config_.unmerged_export_path.empty()) { + if (observer) observer->OnPhase("Writing unmerged reflections"); + WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, + experiment_, true, config_.unmerged_export_path); + logger.Info("Unmerged observations written to {}", config_.unmerged_export_path); + } + if (!config_.unmerged_partials_export_path.empty()) { + if (observer) observer->OnPhase("Writing unmerged reflections"); + WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, + experiment_, false, config_.unmerged_partials_export_path); + logger.Info("Unmerged partials written to {}", config_.unmerged_partials_export_path); + } } // When merging is on, the _process.h5 is skipped because the merged reflections are the wanted diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 23343f65..3c54800c 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -120,9 +120,14 @@ struct ProcessConfig { std::string observation_dump_path; // If set, the integrated observations are written here as an unmerged MTZ (--export-unmerged), - // for aimless / pointless / careless and anything else that scales unmerged data itself. + // for aimless / pointless / careless and anything else that scales unmerged data itself - each + // rocking event summed into one full. std::string unmerged_export_path; + // The same file with each partial written as its own row instead, for the reader to sum + // (--export-unmerged-partials). Independent of the above; a run can write both. + std::string unmerged_partials_export_path; + // If set, after merging, validate the merged intensities against this atomic model (PDB): // compute R-free with an optimized bulk solvent and write 2Fo-Fc / Fo-Fc maps (--model). std::string model_path; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ee8bc783..95f87899 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -162,7 +162,8 @@ void print_usage() { std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; std::cout << " --reference-column