diff --git a/docs/RUGNUX_CLI.md b/docs/RUGNUX_CLI.md index 21d1f1f65..45674058f 100644 --- a/docs/RUGNUX_CLI.md +++ b/docs/RUGNUX_CLI.md @@ -114,8 +114,10 @@ scaled and merged together: a per-frame scale within each sweep, the sweeps' dos their relative B factors (not from the input order or the header times), a decay along that order or one B per sweep, whichever predicts held-out frames better, then the space-group search, the resolution cut and everything after it on the joint data, into `.*`. The detector-dependent -corrections (absorption and detector surfaces) are each sweep's own; the joint merge writes no -unmerged MTZ yet. `--sweeps separate` processes each sweep alone instead, as if it were the only +corrections (absorption and detector surfaces) are each sweep's own. The unmerged MTZ holds the joint +merge's scaled observations in the batches of their own sweeps (sweep `s` of the merge: batches +`s * base + frame + 1`, `base` the power of ten past the longest sweep), under the batch headers that +sweep's own run wrote; the report names each sweep in section `1b. SWEEPS`. `--sweeps separate` processes each sweep alone instead, as if it were the only input, into `_s_*` (k = 1, 2, ... in the order given). What each format needs is in [What Rugnux reads](RUGNUX_FORMATS.md); serial XFEL stills are read diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 79ddbfb6b..5dd08bb6d 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -120,7 +120,8 @@ The header block above section 1 records **how the result was produced**: `RUGNU Rates, per-image costs and progress remain on stdout only. Sections, in order: the `SUMMARY` ([below](#the-summary-and-what-the-run-decided)), then -`1. DATA SET AND GEOMETRY`, `2. CRYSTAL`, `3. MERGED DATA`, `4. DIAGNOSTICS` and, only where a model +`1. DATA SET AND GEOMETRY`, on a joint merge of several sweeps `1b. SWEEPS` +([below](#several-sweeps-merged-as-one)), `2. CRYSTAL`, `3. MERGED DATA`, `4. DIAGNOSTICS` and, only where a model was given (`--model`, `--mr`, or a model file among the inputs), `5. MODEL VALIDATION` - which also holds the keys of [molecular replacement](#molecular-replacement) and of the [refinement](#refinement-refine_) where they ran. A `--mode refine` report has one section, @@ -252,7 +253,7 @@ and a run that produced a textbook data set read identically for their first thr `GONIO_SCALE`, `SPINDLE_CAP`, `ANISOTROPY`, `TWINNING`, `PSEUDO_TRANSLATION`, `LATTICE_TRANSLATION`, `MODEL_HAND`, `MODEL_NOT_VALIDATED`, `MODEL_CELL_MISMATCH`, `MODEL_FREE_SET`, `MODEL_SPACE_GROUP`, `REFERENCE_MISMATCH`, `CANCELLED`, `RESOLUTION_FIT`, `FLIGHT_PATH`, `GEOMETRY_NOT_CONVERGED`, `SCALING_NOT_CONVERGED`, `HARMONIC_CONTAMINATION`, - `SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `LATTICE_EXPLAINS_FEW_SPOTS`, `ICE_RINGS`, `POWDER_RINGS`, `BEAM_CENTRE_OFF`, `OVERLOADS`, `OVERLAPS`. `NONE` when nothing fired. A code appears if and only if its + `SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `LATTICE_EXPLAINS_FEW_SPOTS`, `ICE_RINGS`, `POWDER_RINGS`, `BEAM_CENTRE_OFF`, `OVERLOADS`, `OVERLAPS`, `SWEEP_LEFT_OUT`. `NONE` when nothing fired. A code appears if and only if its warning fired, so the flags and the `WARNING:` lines are two renderings of one list. Thresholds are deliberately low: a warning is a prompt to check something, not a verdict, and some of them fire on data that turn out fine — the closed @@ -286,6 +287,29 @@ and a run that produced a textbook data set read identically for their first thr Then the numbered sections: **1. DATA SET AND GEOMETRY**, **2. CRYSTAL**, **3. MERGED DATA**, **4. DIAGNOSTICS**, and **5. MODEL VALIDATION** only where a model was given. +## Several sweeps merged as one + +A run of several sweeps merged as one (`--sweeps joint`, the default with several rotation inputs, see +[Command-line options](RUGNUX_CLI.md)) describes the joint merge everywhere, and adds section +**1b. SWEEPS**; a run of one input has no such section and no such keys. + +- **`SWEEP_COUNT=`** the sweeps given, **`SWEEP_REFERENCE=`** the one (1-based) the others were put on + the axes of - the one with the most reflections - and **`SWEEP_DECAY_MODEL=`** how the decay over the + sweeps was modelled: `dose-linear` (a B rising along the dose order), `per-sweep B` (one B a sweep) or + `none`, whichever predicted held-out frames better. +- Per sweep `k`: **`SWEEP_k_INPUT=`**, **`SWEEP_k_IMAGES=`**, **`SWEEP_k_INCLUDED=`**; what the sweep + gave processed on its own - **`SWEEP_k_SPACE_GROUP_ALONE=`**, **`SWEEP_k_CELL=`**, + **`SWEEP_k_CC_HALF_ALONE=`**, **`SWEEP_k_COMPLETENESS_ALONE=`**, **`SWEEP_k_RESOLUTION_ALONE=`**; and, + for a sweep in the joint merge, **`SWEEP_k_REINDEX_OP=`** (onto the reference's axes) with + **`SWEEP_k_REINDEX_CC=`** (the merged-intensity correlation under it), **`SWEEP_k_BATCH_FIRST=`** / + **`SWEEP_k_BATCH_LAST=`** (its batches in the unmerged MTZ: sweep `s` of the merge is + `s * base + frame + 1`, `base` the power of ten past the longest sweep), **`SWEEP_k_SCALE=`** (the + median scale of its frames in the joint merge), **`SWEEP_k_RELATIVE_B=`** (its B against the joint + merge, the least damaged sweep at 0) and **`SWEEP_k_DOSE_RANK=`** (its place, from 1, in the dose + order those B factors decide - not the order the sweeps were given in, nor their header times). +- A sweep given but not merged - it gave no merge, or its intensities correlate with the reference + sweep's on no axes - is named by a `SWEEP_LEFT_OUT` warning. + ## The reference-range table `--report-resolution [,]` adds a second block of merging statistics to section 3, the diff --git a/rugnux/ReportDocument.h b/rugnux/ReportDocument.h index 05857375e..2d6d81ecd 100644 --- a/rugnux/ReportDocument.h +++ b/rugnux/ReportDocument.h @@ -27,6 +27,7 @@ namespace PathologyCode { constexpr const char *SYMMETRY_AMBIGUITY = "SYMMETRY_AMBIGUITY"; // the data do not pick a point group constexpr const char *CENTERING_UNTESTED = "CENTERING_UNTESTED"; // centring from the metric only constexpr const char *UNUSABLE_MERGE = "UNUSABLE_MERGE"; // the merged data carry no signal + constexpr const char *SWEEP_LEFT_OUT = "SWEEP_LEFT_OUT"; // a sweep given is not in the joint merge constexpr const char *LOW_COMPLETENESS = "LOW_COMPLETENESS"; // too little of reciprocal space constexpr const char *SWEEP_GAPS = "SWEEP_GAPS"; // ranges that delivered much less constexpr const char *GONIO_SCALE = "GONIO_SCALE"; // the stage rotation is mis-scaled diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 06bc30f0f..57ad0c8bb 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -941,6 +941,64 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " the post-refined geometry. Every number in this report describes the pass named above,\n" " whose files are {}.*; the header-geometry pass is not written.", output_prefix), true)); doc.sections.push_back(std::move(s)); + + } + + // ---------------------------------------------- 1b. SWEEPS (a joint merge of several) + if (!result.joint_sweeps.empty()) { + ReportSection s; + s.title = "1b. SWEEPS"; + Add(s, KeyInt("SWEEP_COUNT", static_cast(result.joint_sweeps.size()))); + Add(s, KeyInt("SWEEP_REFERENCE", result.joint_reference_sweep + 1)); + if (!result.segment_decay_model.empty()) + Add(s, KeyText("SWEEP_DECAY_MODEL", result.segment_decay_model)); + for (size_t k = 0; k < result.joint_sweeps.size(); ++k) { + const auto &w = result.joint_sweeps[k]; + const std::string p = "SWEEP_" + std::to_string(k + 1) + "_"; + const auto key = [&p](const char *name) { return p + name; }; + std::vector names; // the keys outlive the calls that take their c_str() + auto text = [&](const char *name, const std::string &v) { + names.push_back(key(name)); + Add(s, KeyText(names.back().c_str(), v)); + }; + auto real = [&](const char *name, double v, const char *fmt_) { + names.push_back(key(name)); + Add(s, KeyReal(names.back().c_str(), v, fmt_)); + }; + auto integer = [&](const char *name, int64_t v) { + names.push_back(key(name)); + Add(s, KeyInt(names.back().c_str(), v)); + }; + text("INPUT", w.input); + integer("IMAGES", w.images); + text("INCLUDED", w.included ? "TRUE" : "FALSE"); + if (!w.space_group_alone.empty()) + text("SPACE_GROUP_ALONE", w.space_group_alone); + if (w.cell_alone) + text("CELL", CellString(*w.cell_alone)); + if (std::isfinite(w.cc_half_alone)) + real("CC_HALF_ALONE", w.cc_half_alone, "{:.4f}"); + if (std::isfinite(w.completeness_alone)) + real("COMPLETENESS_ALONE", w.completeness_alone, "{:.1f}"); + if (std::isfinite(w.d_min_alone)) + real("RESOLUTION_ALONE", w.d_min_alone, "{:.2f}"); + if (w.included) { + text("REINDEX_OP", w.reindex_op); + real("REINDEX_CC", w.reindex_cc, "{:.4f}"); + integer("BATCH_FIRST", w.batch_first); + integer("BATCH_LAST", w.batch_last); + if (w.segment >= 0 && w.segment < static_cast(result.segment_results.size())) { + const auto &seg = result.segment_results[w.segment]; + if (std::isfinite(seg.scale)) + real("SCALE", seg.scale, "{:.4f}"); + if (std::isfinite(seg.relative_b)) + real("RELATIVE_B", seg.relative_b, "{:.2f}"); + if (seg.dose_rank >= 0) + integer("DOSE_RANK", seg.dose_rank + 1); + } + } + } + doc.sections.push_back(std::move(s)); } // ------------------------------------------------- XDS GEOMETRY KEYWORDS diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 21d77bf6d..223ac20c9 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -913,6 +913,26 @@ struct ProcessResult { // spindle, an indexing ambiguity). The same messages the log warns about, each carrying a closed // PathologyCode for machinery beside the open sentence for a person. std::vector warnings; + // A joint merge of several sweeps (Rugnux::RunJoint, --sweeps joint): what its merge measured per sweep + // segment (RotationScaleMerge::Result::segments, in segment order) and how it modelled the decay, and + // what the run knows of each sweep (filled by the caller, in input order). Empty on a run of one input. + struct JointSegment { double scale = NAN, relative_b = NAN; int dose_rank = -1; }; + std::vector segment_results; + std::string segment_decay_model; + struct JointSweep { + std::string input; + int64_t images = 0; + bool included = false; + int segment = -1; // its segment of the joint merge, where it is in it + std::string space_group_alone; + std::optional cell_alone; + std::string reindex_op; // onto the reference sweep's axes + double reindex_cc = NAN; + int64_t batch_first = 0, batch_last = -1; + double cc_half_alone = NAN, completeness_alone = NAN, d_min_alone = NAN; + }; + std::vector joint_sweeps; + int joint_reference_sweep = -1; // How the crystal sat on the spindle: the angle to the nearest symmetry axis and that axis's // order. A property of the mounting rather than of the data, and the one number that says whether diff --git a/rugnux/RugnuxJointSweeps.cpp b/rugnux/RugnuxJointSweeps.cpp index 09898d4c8..e56ad0066 100644 --- a/rugnux/RugnuxJointSweeps.cpp +++ b/rugnux/RugnuxJointSweeps.cpp @@ -5,9 +5,12 @@ #include #include +#include +#include #include "../common/CrystalLattice.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../gemmi_gph/gemmi/mtz.hpp" JointSweepRelation RelateJointSweeps(const std::vector &sweeps) { JointSweepRelation out; @@ -87,3 +90,93 @@ JointMergeInput BuildJointMergeInput(const std::vector &sweeps } return out; } + +int64_t JointBatchBase(const JointMergeInput &input) { + int64_t longest = 0; + for (const auto &seg : input.segments) + longest = std::max(longest, seg.n_frames); + int64_t base = 1; + while (base <= longest) + base *= 10; + return base; +} + +MtzTable JointUnmergedMtzTable(const std::vector &fulls, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const std::vector &sweeps, + const JointMergeInput &input, float wavelength_A, + const std::vector &history) { + MtzTable mtz; + mtz.space_group_name = space_group.xhm(); + mtz.unit_cell = cell; + mtz.title = "Rugnux unmerged observations, scaled, joint merge of several sweeps"; + mtz.history = history; + mtz.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f}, {"unmerged", "", "unmerged", wavelength_A}}; + for (const auto &[label, type] : std::vector>{ + {"H", 'H'}, {"K", 'H'}, {"L", 'H'}, {"M/ISYM", 'Y'}, {"BATCH", 'B'}, {"I", 'J'}, {"SIGI", 'Q'}, + {"FRACTIONCALC", 'R'}, {"ROT", 'R'}}) + mtz.columns.push_back({label, type, 0, {}}); + + const int64_t base = JointBatchBase(input); + // Per segment, a header for every frame: the sweep's own where its file has one, else its nearest one + // turned by the frames between them (a sweep's file holds the frames its observations span, and a + // batch with no header stops a reading program). + std::vector> own(input.segments.size()); + for (size_t s = 0; s < input.segments.size(); ++s) { + const auto &batches = sweeps[input.segment_sweep[s]].batches; + if (batches.empty()) + continue; + for (int32_t n = 1; n <= input.segments[s].n_frames; ++n) { + const auto nearest = std::min_element(batches.begin(), batches.end(), [n](const auto &a, const auto &b) { + return std::abs(a.number - n) < std::abs(b.number - n); + }); + MtzTable::Batch copy = *nearest; + if (copy.number != n && copy.floats.size() > 47) { + const float shift = static_cast(n - copy.number) * copy.floats[47]; + copy.floats[36] += shift; + copy.floats[37] += shift; + } + own[s][n] = copy; + copy.number = static_cast((s + 1) * base + n); + mtz.batches.push_back(std::move(copy)); + } + } + + const gemmi::GroupOps ops = space_group.operations(); + const auto centring_absent = [&](int h, int k, int l) { + for (const auto &t : ops.cen_ops) + if ((h * t[0] + k * t[1] + l * t[2]) % gemmi::Op::DEN != 0) + return true; + return false; + }; + gemmi::UnmergedHklMover mover(&space_group); + std::vector> rows; + for (const ScaledFull &f : fulls) { + if (centring_absent(f.h, f.k, f.l)) + continue; + const auto frame = static_cast(std::lround(f.image_number)); + size_t s = 0; + while (s + 1 < input.segments.size() && frame >= input.segments[s + 1].first_frame) + ++s; + const int32_t local = frame - input.segments[s].first_frame + 1; // the sweep's own batch + float rot = NAN; + if (const auto it = own[s].find(local); it != own[s].end() && it->second.floats.size() > 47) + rot = it->second.floats[36] + it->second.floats[47] / 2.0f; + std::array hkl{f.h, f.k, f.l}; + const int isym = mover.move_to_asu(hkl); + rows.push_back({static_cast(hkl[0]), static_cast(hkl[1]), static_cast(hkl[2]), + static_cast(isym), static_cast((s + 1) * base + local), f.I, f.sigma, 1.0f, + rot}); + } + std::vector order(rows.size()); + std::iota(order.begin(), order.end(), 0); + std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) { + return std::lexicographical_compare(rows[a].begin(), rows[a].begin() + 5, rows[b].begin(), rows[b].begin() + 5); + }); + for (size_t c = 0; c < mtz.columns.size(); ++c) { + mtz.columns[c].values.resize(rows.size()); + for (size_t i = 0; i < order.size(); ++i) + mtz.columns[c].values[i] = rows[order[i]][c]; + } + mtz.sort_order = {1, 2, 3, 4, 5}; + return mtz; +} diff --git a/rugnux/RugnuxJointSweeps.h b/rugnux/RugnuxJointSweeps.h index d66379830..958d94825 100644 --- a/rugnux/RugnuxJointSweeps.h +++ b/rugnux/RugnuxJointSweeps.h @@ -6,6 +6,7 @@ #include #include +#include "../common/MtzTable.h" #include "../common/Reflection.h" #include "../common/ScalingSettings.h" #include "../common/UnitCell.h" @@ -27,6 +28,8 @@ struct JointMemberSweep { gemmi::SpaceGroup space_group; // the group its run determined std::vector merged; std::vector fulls; // all its fulls, unmerged, on the same axes + std::vector batches; // the batch headers of its own unmerged MTZ + double cc_half = NAN, completeness = NAN, d_min = NAN; // its merge on its own }; // How each sweep is put on the axes of the reference - the sweep with the most fulls: the change of basis @@ -56,3 +59,17 @@ struct JointMergeInput { }; JointMergeInput BuildJointMergeInput(const std::vector &sweeps, const JointSweepRelation &relation); + +// The batch numbers of the joint merge: segment s (from 0) is batches (s + 1) * base + frame + 1, base the +// power of ten past the longest sweep - 1001.. and 2001.. for sweeps of up to 999 frames - so a sweep's +// batches read as its own frames and never run into the next sweep's. +int64_t JointBatchBase(const JointMergeInput &input); + +// The joint merge's unmerged MTZ: its scaled fulls (MergedOutput::scaled_fulls, on the merge's axes and +// frames), each in the batch of its sweep's frame, under the batch headers that sweep's own run wrote - +// its orientation, axis and angles - renumbered (JointBatchBase). Lattice-centring absences are left out, +// as in a single sweep's file. +MtzTable JointUnmergedMtzTable(const std::vector &fulls, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const std::vector &sweeps, + const JointMergeInput &input, float wavelength_A, + const std::vector &history); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 8153d745f..a701fef57 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1671,6 +1671,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.scaling_iterations_partials = r.scaling_iterations_partials; result.scaling_iterations_fulls = r.scaling_iterations_fulls; result.scaling_converged = r.scaling_converged; + result.segment_results.clear(); + for (const auto &seg : r.segments) + result.segment_results.push_back({seg.scale, seg.relative_b, seg.dose_rank}); + result.segment_decay_model = r.segment_decay_model; if (for_search) result.search_scaling_converged = r.scaling_converged; Note(fmt::format("merge in {}: {} unique", label, r.merged.size())); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 32ccf04ae..362539092 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1018,6 +1018,9 @@ struct RugnuxInputRun : RugnuxCli { char centring = 'P'; }; std::optional joint_request; + // What the joint merge adds to its result before anything is written (the sweeps' report keys, the + // unmerged MTZ of its batches). + std::function joint_finish; // Reads the input and sets up the experiment and the run's configuration: RUGNUX_CONTINUE, or the // exit code of a mode that ends here (refine, azint, calibration, scale). @@ -3672,6 +3675,8 @@ int RugnuxInputRun::RunPreparedInput() { throw RugnuxInputFailed{EXIT_FAILURE}; } g_active_process = nullptr; + if (joint_finish) + joint_finish(result); if (sweep_member) { member_result = std::move(result); return EXIT_SUCCESS; @@ -4011,6 +4016,13 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { m.space_group = merged.experiment.GetSpaceGroupOrP1(); m.merged = merged.reflections; m.fulls = merged.scaled_fulls; + if (sweep.member_result->outputs->unmerged) + m.batches = sweep.member_result->outputs->unmerged->batches; + m.cc_half = merged.statistics.overall.cc_half; + m.completeness = merged.statistics.overall.possible_unique_reflections > 0 + ? 100.0 * merged.statistics.overall.unique_reflections + / merged.statistics.overall.possible_unique_reflections : NAN; + m.d_min = merged.statistics.overall.d_min; logger.Info("Sweep {}: {} frames, {} in {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, {} fulls, {} merged", k + 1, m.n_frames, m.space_group.xhm(), m.cell.a, m.cell.b, m.cell.c, m.cell.alpha, m.cell.beta, m.cell.gamma, m.fulls.size(), m.merged.size()); @@ -4063,6 +4075,45 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { ref.space_group.crystal_system(), ref.space_group.centring_type()}; logger.Info("Joint merge of {} sweeps, {} frames", input.segments.size(), joint.joint_request->outcomes.size()); + joint.joint_finish = [&](ProcessResult &result) { + const int64_t base = JointBatchBase(input); + result.joint_reference_sweep = relation.reference; + for (size_t k = 0; k < members.size(); ++k) { + ProcessResult::JointSweep w; + w.input = inputs[k]; + w.images = members[k].n_frames; + w.included = relation.sweeps[k].included; + if (!members[k].merged.empty()) { + w.space_group_alone = members[k].space_group.xhm(); + w.cell_alone = members[k].cell; + w.cc_half_alone = members[k].cc_half; + w.completeness_alone = members[k].completeness; + w.d_min_alone = members[k].d_min; + } + for (size_t s = 0; s < input.segment_sweep.size(); ++s) + if (input.segment_sweep[s] == static_cast(k)) { + w.segment = static_cast(s); + w.batch_first = (s + 1) * base + 1; + w.batch_last = (s + 1) * base + input.segments[s].n_frames; + } + w.reindex_op = relation.sweeps[k].cob.as_hkl().triplet('h'); + w.reindex_cc = relation.sweeps[k].cc; + if (!w.included) + result.warnings.push_back({PathologyCode::SWEEP_LEFT_OUT, fmt::format( + "Sweep {} ({}) is not in the joint merge: {}", k + 1, inputs[k], + relation.sweeps[k].why_left_out.empty() ? std::string("it gave no merge") + : relation.sweeps[k].why_left_out)}); + result.joint_sweeps.push_back(std::move(w)); + } + if (result.outputs && result.outputs->merged && export_unmerged) { + const MergedOutput &merged = *result.outputs->merged; + std::vector history{"From Rugnux " + jfjoch_version() + ", data reduction", + fmt::format("Joint merge of {} sweeps", input.segments.size())}; + result.outputs->unmerged = JointUnmergedMtzTable(merged.scaled_fulls, merged.cell, + merged.experiment.GetSpaceGroupOrP1(), members, input, + merged.experiment.GetWavelength_A(), history); + } + }; rc = joint.RunPreparedInput(); } catch (const RugnuxInputFailed &f) { rc = f.code;