From 5248b1bc009bcd870b5be246119c23ad8b50a231 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 25 Aug 2026 15:21:57 +0200 Subject: [PATCH] rugnux: export the integrated observations as an unmerged MTZ Adds --export-unmerged , in --mode mx and --mode scale alike. rugnux could only write merged reflections, which closed it off from the programs that consume unmerged observations - aimless, pointless, careless and iotbx.merging_statistics all read this file. The column layout follows what pointless itself writes (H K L M/ISYM BATCH I SIGI FRACTIONCALC XDET YDET ROT LP FLAG), with DELPHI, ZETA, BGMEAN and BGVAR added because we have them and a scale model can use them. Three details decide whether the file is usable: Lorentz-polarization is applied and recorded in LP. It is per-observation geometry spanning a factor of ~120 across a sweep and no reading program can recover it; merging with no scale model at all gives R_meas 0.372 / 0.346 / 0.324 in the low shells without it against 0.183 / 0.218 / 0.257 with it. The partiality is not divided out - it stays in FRACTIONCALC - and neither is the per-image scale, since these programs fit their own. The observations are partials and are flagged as such, packed as 256*M + ISYM with LDTYPE=1. Without the flag pointless reads every part as a whole observation and mis-assigns the symmetry; with it, a tetragonal case assembles 8.29M parts into 1.15M observations and comes back as its own space group at confidence 0.92. The scan axis is written as the negation of the stored goniometer axis, which is the Cambridge convention: against pointless's own orientation matrix that agrees to 0.6 degrees, where the axis as stored disagrees by 40. Merging the exported file with aimless reproduces rugnux's own merge to CC 0.9997 on a tetragonal case and 0.9978 on a cubic one, and the anomalous signal survives. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- docs/CHANGELOG.md | 5 + image_analysis/WriteReflections.cpp | 163 ++++++++++++++++++++++++++++ image_analysis/WriteReflections.h | 14 +++ rugnux/Rugnux.cpp | 12 ++ rugnux/Rugnux.h | 4 + rugnux/rugnux_cli.cpp | 16 +++ 6 files changed, 214 insertions(+) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index d2ee03df..b6fe74f3 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,5 +1,10 @@ # Changelog ## 1.0.0 +### 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. + ### 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 20012b5e..cd72a274 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -537,6 +538,168 @@ void WriteShelxHklReflections(const std::vector &reflections, out.close(); } +void WriteUnmergedMtzReflections(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + 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 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(); + 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 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); + } + } + mtz.nreflections = static_cast(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 lattice_at_zero; + std::optional 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(experiment.GetXPixelsNum()); + batch.floats[115] = 1.0f; + batch.floats[116] = static_cast(experiment.GetYPixelsNum()); + + for (const int number : batch_numbers) { + batch.number = number; + batch.floats[36] = phi_start_deg(static_cast(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 &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index fb3f86c5..997b4d97 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -9,6 +9,7 @@ #include "../common/Reflection.h" #include "../common/UnitCell.h" #include "../common/DiffractionExperiment.h" +#include "IntegrationOutcome.h" struct MergeStatistics; struct TwinningAnalysisResult; @@ -46,6 +47,19 @@ 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. +void WriteUnmergedMtzReflections(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename); + void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index fff06f75..00d827c1 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2986,6 +2986,18 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } } + // Unmerged observations (--export-unmerged). After the merge, because that is where the space + // 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); + } + // When merging is on, the _process.h5 is skipped because the merged reflections are the wanted // output and the per-image file is large. If nothing indexed there are no merged reflections, so // the run would finish successfully having written no file at all. Write one now: the per-image diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2e6c40fe..23343f65 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -119,6 +119,10 @@ struct ProcessConfig { // Diagnostic: if set, the -P rot3d combine writes the unmerged fulls here (for comparison vs XDS). 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. + std::string unmerged_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 8217882e..ee8bc783 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -162,6 +162,7 @@ void print_usage() { std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; std::cout << " --reference-column