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