From 4345764a927c04ad0f3aa39600b077129e9f9110 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 12 Nov 2026 06:53:57 +0100 Subject: [PATCH 1/4] rugnux: joint sweeps give the crystal's cell to the sweeps after the first measured one Several sweeps of one crystal of a very large cell (~600 A axis) lost up to half of their sweeps from the joint merge: de novo indexing of a single sweep missed the lattice (an axis harmonic, a P1 supercell, a 5 A row, or no lattice at all), the sweep merged nothing or nothing that correlates with the reference, and was left out. Given the cell the other sweeps measured (-C), the same sweeps index, integrate and merge (checked on five of them, 5 A, ISa 5-10). So once a sweep has merged with a resolved ISa, its cell is given to the sweeps after it as the reference cell; a sweep processed before that which still finds no lattice of the crystal is run again given the shared-cell target's cell. Whether any of them joins is decided on its intensities as before. A user's -C is given to every sweep as before and replaces nothing. Also cheaper relating of the sweeps: the held-out merge of the joined sweeps is built once for all the doubtful sweeps (and once for a re-run sweep's before/after pair) instead of once per sweep, with identical results; the relation is made a second time only where a sweep's run was replaced; TIMING marks for its stages. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- rugnux/RugnuxJointSweeps.cpp | 204 +++++++++++++++++++++-------------- rugnux/RugnuxJointSweeps.h | 15 ++- rugnux/rugnux_cli.cpp | 49 +++++++-- 3 files changed, 175 insertions(+), 93 deletions(-) diff --git a/rugnux/RugnuxJointSweeps.cpp b/rugnux/RugnuxJointSweeps.cpp index 6da5805fd..221893b1e 100644 --- a/rugnux/RugnuxJointSweeps.cpp +++ b/rugnux/RugnuxJointSweeps.cpp @@ -11,6 +11,7 @@ #include #include "../common/CrystalLattice.h" +#include "../common/Logger.h" #include #include "../image_analysis/scale_merge/HKLKey.h" @@ -109,10 +110,11 @@ SweepIsomorphism MeasureSweepIsomorphism(const std::vector &a, const // The merge with and without a data set, after the delta-CC1/2 of Assmann, Brehm & Diederichs (2016) // J. Appl. Cryst. 49, 1021-1028 -HeldOutGain MeasureHeldOutGain(const std::vector *> &base, - const std::vector &base_cob, - const std::vector &d, const gemmi::Op &cob_d, - const gemmi::SpaceGroup &space_group) { +std::vector MeasureHeldOutGains(const std::vector *> &base, + const std::vector &base_cob, + const std::vector *> &ds, + const std::vector &cob_ds, + const gemmi::SpaceGroup &space_group) { const HKLKeyGenerator key(/*merge_friedel=*/true, space_group); struct Obs { float I, sigma, d; }; const auto group = [&](const std::vector &fulls, const gemmi::Op &cob) { @@ -180,84 +182,98 @@ HeldOutGain MeasureHeldOutGain(const std::vector * for (auto &[k, obs] : m) pool[k].insert(pool[k].end(), obs.begin(), obs.end()); } - auto added = group(d, cob_d); - put_on_scale(added, reference); + // The merge of `base` is the same for every sweep measured against it, so it is made once. + const auto measure = [&](const std::vector &d, const gemmi::Op &cob_d) { + auto added = group(d, cob_d); + put_on_scale(added, reference); - // Per reflection both have, with observations of the merge in both halves: the weighted squared - // residual of each half against the other half's mean, without and with the doubtful sweep's. - struct Row { double without, with; }; // per held-out observation - std::vector rows; - for (const auto &[k, obs] : pool) { - const auto it = added.find(k); - if (it == added.end() || obs.size() < 2) - continue; - double wd = 0, wid = 0; - for (const Obs &o : it->second) { - wd += 1.0 / (static_cast(o.sigma) * o.sigma); - wid += o.I / (static_cast(o.sigma) * o.sigma); - } - Row r{0, 0}; - double n_target = 0; - for (int h = 0; h < 2; ++h) { - double w = 0, wi = 0; // the predicting half - for (size_t i = 1 - h; i < obs.size(); i += 2) { - w += 1.0 / (static_cast(obs[i].sigma) * obs[i].sigma); - wi += obs[i].I / (static_cast(obs[i].sigma) * obs[i].sigma); + // Per reflection both have, with observations of the merge in both halves: the weighted squared + // residual of each half against the other half's mean, without and with the doubtful sweep's. + struct Row { double without, with; }; // per held-out observation + std::vector rows; + for (const auto &[k, obs] : pool) { + const auto it = added.find(k); + if (it == added.end() || obs.size() < 2) + continue; + double wd = 0, wid = 0; + for (const Obs &o : it->second) { + wd += 1.0 / (static_cast(o.sigma) * o.sigma); + wid += o.I / (static_cast(o.sigma) * o.sigma); } - const double without = wi / w, with = (wi + wid) / (w + wd); - for (size_t i = h; i < obs.size(); i += 2) { - const double wt = 1.0 / (static_cast(obs[i].sigma) * obs[i].sigma); - r.without += wt * (obs[i].I - without) * (obs[i].I - without); - r.with += wt * (obs[i].I - with) * (obs[i].I - with); - n_target += 1; + Row r{0, 0}; + double n_target = 0; + for (int h = 0; h < 2; ++h) { + double w = 0, wi = 0; // the predicting half + for (size_t i = 1 - h; i < obs.size(); i += 2) { + w += 1.0 / (static_cast(obs[i].sigma) * obs[i].sigma); + wi += obs[i].I / (static_cast(obs[i].sigma) * obs[i].sigma); + } + const double without = wi / w, with = (wi + wid) / (w + wd); + for (size_t i = h; i < obs.size(); i += 2) { + const double wt = 1.0 / (static_cast(obs[i].sigma) * obs[i].sigma); + r.without += wt * (obs[i].I - without) * (obs[i].I - without); + r.with += wt * (obs[i].I - with) * (obs[i].I - with); + n_target += 1; + } } + rows.push_back({r.without / n_target, r.with / n_target}); } - rows.push_back({r.without / n_target, r.with / n_target}); - } - HeldOutGain out; - out.base_unique = static_cast(pool.size()); - for (const auto &[k, obs] : added) - if (!pool.contains(k)) - ++out.added; - constexpr size_t MIN_ROWS = 50; - if (rows.size() < MIN_ROWS) - return out; - // Medians, not sums: one reflection a twin overlap or an ice ring spoils would otherwise decide it. - // The median change per reflection, relative to the median residual without the sweep. - const auto change = [&](const std::vector &idx) { - std::vector diff(idx.size()), without(idx.size()); - for (size_t j = 0; j < idx.size(); ++j) { - diff[j] = rows[idx[j]].with - rows[idx[j]].without; - without[j] = rows[idx[j]].without; - } - const auto median = [](std::vector &v) { - std::nth_element(v.begin(), v.begin() + v.size() / 2, v.end()); - return v[v.size() / 2]; + HeldOutGain out; + out.base_unique = static_cast(pool.size()); + for (const auto &[k, obs] : added) + if (!pool.contains(k)) + ++out.added; + constexpr size_t MIN_ROWS = 50; + if (rows.size() < MIN_ROWS) + return out; + // Medians, not sums: one reflection a twin overlap or an ice ring spoils would otherwise decide it. + // The median change per reflection, relative to the median residual without the sweep. + const auto change = [&](const std::vector &idx) { + std::vector diff(idx.size()), without(idx.size()); + for (size_t j = 0; j < idx.size(); ++j) { + diff[j] = rows[idx[j]].with - rows[idx[j]].without; + without[j] = rows[idx[j]].without; + } + const auto median = [](std::vector &v) { + std::nth_element(v.begin(), v.begin() + v.size() / 2, v.end()); + return v[v.size() / 2]; + }; + const double base = median(without); + return base > 0 ? median(diff) / base : NAN; }; - const double base = median(without); - return base > 0 ? median(diff) / base : NAN; + std::vector idx(rows.size()); + std::iota(idx.begin(), idx.end(), 0); + out.change = change(idx); + out.n = static_cast(rows.size()); + constexpr int BOOTSTRAP = 200; + std::mt19937 rng(12345); + std::uniform_int_distribution pick(0, rows.size() - 1); + double s = 0, ss = 0; + int m = 0; + for (int t = 0; t < BOOTSTRAP; ++t) { + for (size_t &i : idx) + i = pick(rng); + const double c = change(idx); + if (!std::isfinite(c)) continue; + s += c; ss += c * c; ++m; + } + if (m > 1) + out.se = std::sqrt(std::max(0.0, (ss - s * s / m) / (m - 1))); + return out; }; - std::vector idx(rows.size()); - std::iota(idx.begin(), idx.end(), 0); - out.change = change(idx); - out.n = static_cast(rows.size()); - constexpr int BOOTSTRAP = 200; - std::mt19937 rng(12345); - std::uniform_int_distribution pick(0, rows.size() - 1); - double s = 0, ss = 0; - int m = 0; - for (int t = 0; t < BOOTSTRAP; ++t) { - for (size_t &i : idx) - i = pick(rng); - const double c = change(idx); - if (!std::isfinite(c)) continue; - s += c; ss += c * c; ++m; - } - if (m > 1) - out.se = std::sqrt(std::max(0.0, (ss - s * s / m) / (m - 1))); + std::vector out; + for (size_t i = 0; i < ds.size(); ++i) + out.push_back(measure(*ds[i], cob_ds[i])); return out; } +HeldOutGain MeasureHeldOutGain(const std::vector *> &base, + const std::vector &base_cob, + const std::vector &d, const gemmi::Op &cob_d, + const gemmi::SpaceGroup &space_group) { + return MeasureHeldOutGains(base, base_cob, {&d}, {cob_d}, space_group)[0]; +} + namespace { // h -> (h cob) law, as one operator on Miller indices (apply_to_hkl is h R / DEN). gemmi::Op ThenLaw(const gemmi::Op &cob, const gemmi::Op &law) { @@ -390,6 +406,8 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps // First the axes: each sweep on the reference's (ChooseReferenceSetting), then the twin laws of the // lattice settled over all the pairs at once (ConsensusTwinLaws). + Logger logger("Rugnux"); + TimingMark axes_mark(logger, "joint sweeps: axes"); for (size_t k = 0; k < sweeps.size(); ++k) { if (static_cast(k) == out.reference) continue; @@ -410,8 +428,12 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps s.why_left_out = "it shares too few reflections with the reference sweep to be put on its axes"; } } + axes_mark.Stop(); + TimingMark laws_mark(logger, "joint sweeps: twin laws"); ConsensusTwinLaws(sweeps, out); + laws_mark.Stop(); + TimingMark iso_mark(logger, "joint sweeps: isomorphism"); std::vector doubtful; for (size_t k = 0; k < sweeps.size(); ++k) { JointSweepRelation::Sweep &s = out.sweeps[k]; @@ -447,7 +469,10 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps doubtful.push_back(k); } + iso_mark.Stop(); + // The doubtful sweeps, each against the merge of the clearly joined ones. + TimingMark held_out_mark(logger, "joint sweeps: held-out"); std::vector *> base{&ref.fulls}; std::vector base_cob{gemmi::Op::identity()}; for (size_t k = 0; k < sweeps.size(); ++k) @@ -455,9 +480,18 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps base.push_back(&sweeps[k].fulls); base_cob.push_back(out.sweeps[k].cob); } + std::vector *> doubtful_fulls; + std::vector doubtful_cob; for (const size_t k : doubtful) { + doubtful_fulls.push_back(&sweeps[k].fulls); + doubtful_cob.push_back(out.sweeps[k].cob); + } + const std::vector gains = MeasureHeldOutGains(base, base_cob, doubtful_fulls, doubtful_cob, + ref.space_group); + for (size_t i = 0; i < doubtful.size(); ++i) { + const size_t k = doubtful[i]; JointSweepRelation::Sweep &s = out.sweeps[k]; - const HeldOutGain g = MeasureHeldOutGain(base, base_cob, sweeps[k].fulls, s.cob, ref.space_group); + const HeldOutGain &g = gains[i]; s.held_out_change = g.change; s.held_out_se = g.se; s.held_out_n = g.n; @@ -478,16 +512,24 @@ double VolumeLength(const UnitCell &cell) { return std::cbrt(static_cast(CrystalLattice(cell).CalcVolume())); } -HeldOutGain HeldOutAgainst(const std::vector &base, const std::vector &base_cob, - const JointMemberSweep &reference, const JointMemberSweep &d) { - const ReferenceSettingChoice choice = ChooseReferenceSetting(d.merged, d.cell, d.space_group, reference.merged, - reference.cell, reference.space_group); - if (!std::isfinite(choice.cc)) - return {}; +std::vector HeldOutAgainst(const std::vector &base, + const std::vector &base_cob, const JointMemberSweep &reference, + const std::vector &ds) { std::vector *> fulls; for (const JointMemberSweep *m : base) fulls.push_back(&m->fulls); - return MeasureHeldOutGain(fulls, base_cob, d.fulls, choice.cob, reference.space_group); + // Each on the reference's axes; one that cannot be put there is measured as having nothing in common. + static const std::vector none; + std::vector *> d_fulls; + std::vector d_cob; + for (const JointMemberSweep *d : ds) { + const ReferenceSettingChoice choice = ChooseReferenceSetting(d->merged, d->cell, d->space_group, + reference.merged, reference.cell, + reference.space_group); + d_fulls.push_back(std::isfinite(choice.cc) ? &d->fulls : &none); + d_cob.push_back(choice.cob); + } + return MeasureHeldOutGains(fulls, base_cob, d_fulls, d_cob, reference.space_group); } int SharedCellTarget(const std::vector &sweeps, const JointSweepRelation &relation) { diff --git a/rugnux/RugnuxJointSweeps.h b/rugnux/RugnuxJointSweeps.h index 4150100c9..dbc0afd8c 100644 --- a/rugnux/RugnuxJointSweeps.h +++ b/rugnux/RugnuxJointSweeps.h @@ -75,6 +75,12 @@ HeldOutGain MeasureHeldOutGain(const std::vector * const std::vector &base_cob, const std::vector &d, const gemmi::Op &cob_d, const gemmi::SpaceGroup &space_group); +// The same for several sweeps `ds`, each alone against the one merge of `base`. +std::vector MeasureHeldOutGains(const std::vector *> &base, + const std::vector &base_cob, + const std::vector *> &ds, + const std::vector &cob_ds, + const gemmi::SpaceGroup &space_group); // How each sweep is put on the axes of the reference and whether it joins it. The reference is the // best-diffracting sweep whose intensities have a measured accuracy: of the sweeps whose error model @@ -150,10 +156,11 @@ void ConsensusTwinLaws(const std::vector &sweeps, JointSweepRe // The cube root of the cell's volume: a length to compare the cells of sweeps on different axes by. double VolumeLength(const UnitCell &cell); -// MeasureHeldOutGain of sweep `d` against the merge of `base` (on `reference`'s axes by `base_cob`, the -// first of them the scale), `d` put on `reference`'s axes first (ChooseReferenceSetting). -HeldOutGain HeldOutAgainst(const std::vector &base, const std::vector &base_cob, - const JointMemberSweep &reference, const JointMemberSweep &d); +// MeasureHeldOutGains of the sweeps `ds` against the merge of `base` (on `reference`'s axes by `base_cob`, the +// first of them the scale), each put on `reference`'s axes first (ChooseReferenceSetting). +std::vector HeldOutAgainst(const std::vector &base, + const std::vector &base_cob, const JointMemberSweep &reference, + const std::vector &ds); // The sweep whose cell the crystal's is taken to be: of the sweeps that join, the one reaching the finest // resolution on its own (ties: the higher ISa). A cell and the detector distance are one degeneracy at low diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index fe6dea69c..87cc2278b 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -4094,13 +4094,16 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { } } // Each sweep processed on its own first, writing nothing, its fulls kept: the first level. `distance_mm` - // replaces the file's detector distance where given. - const auto run_member = [&](size_t k, std::optional distance_mm, const std::string &log_suffix) { + // replaces the file's detector distance where given, and `cell` is the reference cell (-C) where given. + const auto run_member = [&](size_t k, std::optional distance_mm, std::optional cell, + const std::string &log_suffix) { RugnuxInputRun sweep(*this); sweep.input_files = {plan.sweeps[k].input}; sweep.log_open = false; if (distance_mm) sweep.detector_distance_mm = distance_mm; + if (cell) + sweep.fixed_reference_unit_cell = cell; // A sweep of a joint run writes its log, and nothing else (sweep_member). if (!output_prefix.empty()) sweep.output_prefix = output_prefix + "_s" + std::to_string(k + 1) + log_suffix; @@ -4149,18 +4152,46 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { }; std::vector members; std::vector inputs; + // The sweeps are of one crystal: once one has merged with a measured accuracy (its error model resolved + // an ISa, as the reference is chosen), its cell is given to the sweeps after it as a reference cell (-C), + // which de novo indexing of one sweep can otherwise miss on a large cell - and a lattice it misses costs + // more than a whole run before it is found to merge nothing. Where -C was given, every sweep has it. + std::optional crystal_cell; + std::vector cell_given(plan.sweeps.size(), false); for (size_t k = 0; k < plan.sweeps.size(); k++) { logger.Info("Sweep {} of {}: {}, processed on its own", k + 1, plan.sweeps.size(), plan.sweeps[k].label); - members.push_back(run_member(k, std::nullopt, "")); + cell_given[k] = crystal_cell.has_value() || fixed_reference_unit_cell.has_value(); + members.push_back(run_member(k, std::nullopt, crystal_cell, "")); inputs.push_back(plan.sweeps[k].input); + if (!cell_given[k] && members[k].isa > 0.0) { + crystal_cell = members[k].cell; + logger.Info("Sweep {}: its cell is given to the sweeps after it", k + 1); + } } // One cell for the crystal (SharedCellTarget): a sweep whose own cell disagrees with it - joined or not // on the first look - is processed again from the distance that would make the two agree, and that run replaces its first one where it // agrees better with the target sweep's intensities (lower held-out residual, HeldOutAgainst). + JointSweepRelation relation = RelateJointSweeps(members); + bool replaced = false; // a sweep's run replaced by its run again: the relation is made again { - const JointSweepRelation first = RelateJointSweeps(members); + const JointSweepRelation &first = relation; const int target = SharedCellTarget(members, first); + // A sweep processed before any cell was given that found no lattice of this crystal - it merged + // nothing, or nothing that correlates with the reference on any axes - is processed again given the + // target sweep's cell. Its join is then decided on its intensities, as every sweep's is. + for (size_t k = 0; target >= 0 && k < members.size(); ++k) { + const std::string &decision = first.sweeps[k].decision; + if (cell_given[k] || (decision != "no merge" && decision != "no axes")) + continue; + logger.Info("Sweep {}: {}; processed again given sweep {}'s cell", k + 1, first.sweeps[k].why_left_out, + target + 1); + JointMemberSweep again = run_member(k, std::nullopt, members[target].cell, "_known"); + if (!again.merged.empty()) { + members[k] = std::move(again); + replaced = true; + } + } for (size_t k = 0; target >= 0 && k < members.size(); ++k) { // Every sweep on the reference's axes, joined or not: a cell off by a few percent is itself a reason // for intensities to disagree, so the join is decided after this, on the runs it keeps. @@ -4174,7 +4205,7 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { const auto distance = static_cast(members[k].distance_mm * scale); logger.Info("Sweep {}: its cell is {:+.2f}% off sweep {}'s; processed again from {:.2f} mm", k + 1, 100.0 * (1.0 / scale - 1.0), target + 1, distance); - JointMemberSweep again = run_member(k, distance, "_cell"); + JointMemberSweep again = run_member(k, distance, std::nullopt, "_cell"); if (again.merged.empty()) continue; // Against the merge of the other joined sweeps, the target first (the scale). @@ -4186,8 +4217,8 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { base_cob.push_back(first.sweeps[j].cob); } const JointMemberSweep &reference = members[first.reference]; - const HeldOutGain before = HeldOutAgainst(base, base_cob, reference, members[k]); - const HeldOutGain after = HeldOutAgainst(base, base_cob, reference, again); + const std::vector gains = HeldOutAgainst(base, base_cob, reference, {&members[k], &again}); + const HeldOutGain &before = gains[0], &after = gains[1]; const bool adopt = std::isfinite(after.change) && !(after.change >= before.change); logger.Info("Sweep {}: against the joined sweeps (scale: sweep {}), held-out residual {:+.2f}% at its own cell ({:.3f} {:.3f} {:.3f}), " "{:+.2f}% processed again ({:.3f} {:.3f} {:.3f}): {}", k + 1, target + 1, @@ -4197,12 +4228,14 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { if (adopt) { again.reprocessed_from_mm = members[k].distance_mm; members[k] = std::move(again); + replaced = true; } } } // On one set of axes, and laid end to end: the second level's input. - const JointSweepRelation relation = RelateJointSweeps(members); + if (replaced) + relation = RelateJointSweeps(members); if (relation.reference < 0 || members[relation.reference].fulls.empty()) { logger.Error("No sweep merged on its own, so there is nothing to merge jointly"); return EXIT_FAILURE; From ec3f2ebb800e7e86d2fd00fcab5886407e59b857 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 12 Nov 2026 07:53:25 +0100 Subject: [PATCH 2/4] rugnux: joint sweeps - only a sweep that failed is given the crystal's cell; relations measured in parallel The previous commit gave the first measured sweep's cell to every sweep after it. On a crystal of a ~600 A cell with a pseudo-centring that moved sweeps which index well de novo into another setting or centring (two of twenty lost, two gained: still 13 joined), so the cell now goes only to a sweep that found no lattice of the crystal on its own - it merged nothing, or nothing that correlates with the reference on any axes - and that sweep is processed once more given the shared-cell target's cell. A space group given as well (-S) was tried and refused by the metric check on such sweeps. The isomorphism ratios and the held-out gains of the sweeps are independent of each other and are now measured at once, one task per sweep (the held-out merge of the joined sweeps still built once). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- rugnux/RugnuxJointSweeps.cpp | 24 +++++++++++++++++++----- rugnux/rugnux_cli.cpp | 26 +++++++++----------------- 2 files changed, 28 insertions(+), 22 deletions(-) diff --git a/rugnux/RugnuxJointSweeps.cpp b/rugnux/RugnuxJointSweeps.cpp index 221893b1e..c6a6253e0 100644 --- a/rugnux/RugnuxJointSweeps.cpp +++ b/rugnux/RugnuxJointSweeps.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -261,9 +262,13 @@ std::vector MeasureHeldOutGains(const std::vector out; + // ... and each sweep against it on its own, all of them at once. + std::vector> gains; for (size_t i = 0; i < ds.size(); ++i) - out.push_back(measure(*ds[i], cob_ds[i])); + gains.push_back(std::async(std::launch::async, [&measure, &ds, &cob_ds, i] { return measure(*ds[i], cob_ds[i]); })); + std::vector out; + for (auto &g : gains) + out.push_back(g.get()); return out; } @@ -434,7 +439,8 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps laws_mark.Stop(); TimingMark iso_mark(logger, "joint sweeps: isomorphism"); - std::vector doubtful; + // Each sweep's ratio is measured on its own, all of them at once. + std::vector> isomorphism(sweeps.size()); for (size_t k = 0; k < sweeps.size(); ++k) { JointSweepRelation::Sweep &s = out.sweeps[k]; if (static_cast(k) == out.reference || !s.decision.empty()) @@ -445,8 +451,16 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps "(best CC {:.3f})", s.cc); continue; } - const SweepIsomorphism iso = MeasureSweepIsomorphism(sweeps[k].fulls, s.cob, ref.fulls, - gemmi::Op::identity(), ref.space_group); + isomorphism[k] = std::async(std::launch::async, [&sweeps, &ref, cob = s.cob, k] { + return MeasureSweepIsomorphism(sweeps[k].fulls, cob, ref.fulls, gemmi::Op::identity(), ref.space_group); + }); + } + std::vector doubtful; + for (size_t k = 0; k < sweeps.size(); ++k) { + if (!isomorphism[k].valid()) + continue; + JointSweepRelation::Sweep &s = out.sweeps[k]; + const SweepIsomorphism iso = isomorphism[k].get(); s.iso_ratio = iso.ratio; s.iso_se = iso.se; s.iso_n = iso.n; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 87cc2278b..13589b719 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -4152,21 +4152,10 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { }; std::vector members; std::vector inputs; - // The sweeps are of one crystal: once one has merged with a measured accuracy (its error model resolved - // an ISa, as the reference is chosen), its cell is given to the sweeps after it as a reference cell (-C), - // which de novo indexing of one sweep can otherwise miss on a large cell - and a lattice it misses costs - // more than a whole run before it is found to merge nothing. Where -C was given, every sweep has it. - std::optional crystal_cell; - std::vector cell_given(plan.sweeps.size(), false); for (size_t k = 0; k < plan.sweeps.size(); k++) { logger.Info("Sweep {} of {}: {}, processed on its own", k + 1, plan.sweeps.size(), plan.sweeps[k].label); - cell_given[k] = crystal_cell.has_value() || fixed_reference_unit_cell.has_value(); - members.push_back(run_member(k, std::nullopt, crystal_cell, "")); + members.push_back(run_member(k, std::nullopt, std::nullopt, "")); inputs.push_back(plan.sweeps[k].input); - if (!cell_given[k] && members[k].isa > 0.0) { - crystal_cell = members[k].cell; - logger.Info("Sweep {}: its cell is given to the sweeps after it", k + 1); - } } // One cell for the crystal (SharedCellTarget): a sweep whose own cell disagrees with it - joined or not @@ -4177,12 +4166,15 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { { const JointSweepRelation &first = relation; const int target = SharedCellTarget(members, first); - // A sweep processed before any cell was given that found no lattice of this crystal - it merged - // nothing, or nothing that correlates with the reference on any axes - is processed again given the - // target sweep's cell. Its join is then decided on its intensities, as every sweep's is. - for (size_t k = 0; target >= 0 && k < members.size(); ++k) { + // A sweep that found no lattice of this crystal on its own - it merged nothing, or nothing that + // correlates with the reference on any axes - is processed again given the target sweep's cell (-C): + // de novo indexing of one sweep can miss a large cell that the crystal's other sweeps have measured. + // Its join is then decided on its intensities, as every sweep's is. Only a sweep that failed is given + // the cell: given to every sweep from the start, it moved sweeps that index well de novo into another + // setting or centring. Where -C was given, every sweep already had a cell. + for (size_t k = 0; target >= 0 && !fixed_reference_unit_cell && k < members.size(); ++k) { const std::string &decision = first.sweeps[k].decision; - if (cell_given[k] || (decision != "no merge" && decision != "no axes")) + if (decision != "no merge" && decision != "no axes") continue; logger.Info("Sweep {}: {}; processed again given sweep {}'s cell", k + 1, first.sweeps[k].why_left_out, target + 1); From 59ab1935e4fc9def82e5170fe316f9891b5aac3c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 12 Nov 2026 11:57:14 +0100 Subject: [PATCH 3/4] rugnux: tNCS search bounded on very large merges; joint sweeps: one search arm, no report-only tNCS, refusals - tNCS: the modulation search reads at most the 200,000 coarsest reflections of its band, and the ladder that refines the class vector climbs to at most the 1,000,000 lowest orders. Both grew with the merge: on a ~600 A cell merged in P1 (2-3M reflections) the pre-promotion analysis took 38-46 s per sweep, 23-26 s of it on the critical path; now ~3 s + the ladder on 1M. The tNCS truth sets read up to ~110,000 and ~960,000 and are unchanged (9ea5, 8c3e, 6yqf, 5ebi, 2wnq, 5j23, 6pxb, 9i80, four cytc sets: identical rows); a 6 M-reflection open set (4nwv) keeps its result, 131 -> 96 s. - A sweep of a joint run skips the final tNCS analysis, which is reporting only and which it does not report. - The joint merge searches the space group on one arm: its fulls come in pre-combined at zeta 1, so the arm of only the well-measured observations read all of them again (one P1 merge, ~90 s on 20 sweeps of a large cell, and its engine's memory). - A sweep of a joint run that ends on a JFJochException (its lattice merged nothing, say) is logged as that sweep leaving the joint merge, a warning, not "Processing failed": the joint run goes on and reports it. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- image_analysis/scale_merge/TranslationalNCS.cpp | 17 ++++++++++++++++- rugnux/RugnuxScaleMerge.cpp | 5 +++-- rugnux/rugnux_cli.cpp | 12 ++++++++++++ 3 files changed, 31 insertions(+), 3 deletions(-) diff --git a/image_analysis/scale_merge/TranslationalNCS.cpp b/image_analysis/scale_merge/TranslationalNCS.cpp index a5ccddd5c..c72909778 100644 --- a/image_analysis/scale_merge/TranslationalNCS.cpp +++ b/image_analysis/scale_merge/TranslationalNCS.cpp @@ -46,6 +46,15 @@ namespace { // 4 A, as the Patterson band is: a small cell's |h| at 2.5 A is a protein cell's at 4 A, so // the phase is no more sensitive to an error in u there. constexpr int kMinModulationBand = 2000; + // ... and bounded to the coarsest this many: past it the search reads no more than the precision it + // needs, and its cost, a cosine per reflection per trial vector, grew with the merge (45 s on a + // ~600 A cell merged in P1, 2M reflections to 4 A). The tNCS truth sets read up to ~110,000. + constexpr int kMaxModulationBand = 200000; + // The ladder that refines the class vector on the whole data (below) climbs to at most this many + // reflections, in order of |h|: on a merge of millions (a ~600 A cell in P1) the climb past them took + // 20 s and refined u at orders a merge of that size reaches only at its edge. The tNCS truth sets + // read up to ~960,000. + constexpr size_t kMaxClassReflections = 1000000; constexpr int kModulationShells = 10; constexpr int kPhaseBins = 12; constexpr int kMinModulationReflections = 500; @@ -484,7 +493,8 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector merged[b].d : a < b; }); const int n_beyond = static_cast(std::count_if(mod_idx.begin(), mod_idx.end(), [&](int i) { return merged[i].d >= kModulationDMin; })); - mod_idx.resize(std::max(n_beyond, std::min(kMinModulationBand, static_cast(mod_idx.size())))); + mod_idx.resize(std::min(kMaxModulationBand, std::max(n_beyond, std::min(kMinModulationBand, + static_cast(mod_idx.size()))))); if (!mod_idx.empty()) { // Shells linear in 1/d^2, and the same defensive normalisation the twinning statistics // use: noise-only shells skipped, Wilson outliers rejected, the shell mean re-fitted once @@ -832,6 +842,11 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector kMaxClassReflections) { + cls_hkl.resize(kMaxClassReflections); + cls_e2.resize(kMaxClassReflections); + } } // The steps of the ladder double the number of reflections read rather than the |h| ceiling, // so that the first step always has a population to refine on whatever the cell is and no diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 8b74ffddc..2e270683a 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -4480,9 +4480,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (!geometry_prepass && !superseded) { // Translational pseudo-symmetry first: its verdict is an input to the twinning analysis, // which must not read its own L-test when a non-half-integer pseudo-translation has biased - // it. Reporting only - no reflection and no merge depends on the outcome. + // it. Reporting only - no reflection and no merge depends on the outcome, so a sweep of a joint + // run, which reports nothing of its own, does not measure it. TimingMark tncs_mark(logger, "analyses/tNCS"); - if (result.consensus_cell) + if (result.consensus_cell && !config_.sweep_member) result.tncs = AnalyzeTranslationalNCS(sm.merged, gemmi::UnitCell(*result.consensus_cell), twin_sg, 30, 0, config_.nthreads); tncs_mark.Stop(); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 13589b719..675567094 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -3773,6 +3773,15 @@ int RugnuxInputRun::RunPreparedInput() { console->Fail(); PrintSummary(*console, doc, output_prefix, files_since, provenance.wall_time_s); throw RugnuxInputFailed{EXIT_FAILURE}; + } catch (const JFJochException &e) { + // A sweep of a joint run that ends without a merge (its lattice merged nothing, say) is left out of + // the joint merge, which goes on: the run reports it there, as a refusal of that sweep. + if (sweep_member) + logger.Warning("This sweep stops here and is left out of the joint merge: {}", FailureText(e)); + else + logger.Error("Processing failed: {}", FailureText(e)); + sink.reset(); + throw RugnuxInputFailed{EXIT_FAILURE}; } catch (const std::exception &e) { logger.Error("Processing failed: {}", FailureText(e)); sink.reset(); @@ -4273,6 +4282,9 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { return rc; ScalingSettings s = joint.experiment.GetScalingSettings(); s.SweepSegments(input.segments); + // The fulls come in pre-combined, at zeta 1 (BuildJointMergeInput), so a search arm of the + // well-measured observations only would read all of them again: one search. + s.SearchMinZeta(0.0); joint.experiment.ImportScalingSettings(s); joint.config.rotation_postrefine_geometry = false; joint.config.export_unmerged = false; From 30770ddb1dc2f412b7dd8fa1d307388bd47fcc68 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 12 Nov 2026 11:59:10 +0100 Subject: [PATCH 4/4] rugnux: joint sweeps release the members' fulls, and the joint input holds no spare capacity The members' fulls were "released" with `= {}`, an assignment from an empty initializer list that keeps the vectors' memory (about 0.6 GB per sweep of a ~600 A cell, ~11 GB for 20 sweeps, held through the joint merge); they are now assigned empty vectors. The joint merge's input grew each frame's reflections one by one, up to twice their size in spare capacity; each frame is now reserved at its count first. Nothing the merge reads changes. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- rugnux/RugnuxJointSweeps.cpp | 19 ++++++++++++++++--- rugnux/rugnux_cli.cpp | 5 +++-- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/rugnux/RugnuxJointSweeps.cpp b/rugnux/RugnuxJointSweeps.cpp index c6a6253e0..db0a59943 100644 --- a/rugnux/RugnuxJointSweeps.cpp +++ b/rugnux/RugnuxJointSweeps.cpp @@ -574,11 +574,24 @@ JointMergeInput BuildJointMergeInput(const std::vector &sweeps out.outcomes.resize(first + sweep.n_frames); for (int64_t f = 0; f < sweep.n_frames; ++f) out.outcomes[first + f].latt = lattice; + const auto usable = [](const ScaledFull &full) { + return std::isfinite(full.I) && std::isfinite(full.var_bkg) && full.var_per_I > 0.0f && full.d > 0.0f; + }; + const auto frame_of = [&](const ScaledFull &full) { + return std::clamp(std::llround(full.image_number), 0, sweep.n_frames - 1); + }; + // Each frame's reflections reserved at their count first: grown one by one, a large cell's + // sweeps held up to twice their size in spare capacity. + std::vector count(sweep.n_frames, 0); + for (const ScaledFull &full : sweep.fulls) + if (usable(full)) + count[frame_of(full)]++; + for (int64_t f = 0; f < sweep.n_frames; ++f) + out.outcomes[first + f].reflections.reserve(count[f]); for (const ScaledFull &full : sweep.fulls) { - if (!std::isfinite(full.I) || !std::isfinite(full.var_bkg) || !(full.var_per_I > 0.0f) - || !(full.d > 0.0f)) + if (!usable(full)) continue; - const auto frame = std::clamp(std::llround(full.image_number), 0, sweep.n_frames - 1); + const auto frame = frame_of(full); const gemmi::Op::Miller hkl = rel.cob.apply_to_hkl({{full.h, full.k, full.l}}); // A full as the merge takes a pre-combined observation, so that the combine rebuilds the variance // its own sweep's merge gave it: a count c of intensity I/c at correction c = var_per_I is a diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 675567094..9ca68f76a 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -4262,9 +4262,10 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { JointMergeInput input = BuildJointMergeInput(members, relation); // The joint merge works on its own copy of the fulls (input.outcomes): the members' go, so that a large // cell's sweeps are not held twice. + // (Assigned an empty vector, not {}: that is an assignment from an empty list, which keeps the memory.) for (JointMemberSweep &m : members) { - m.fulls = {}; - m.merged = {}; + m.fulls = std::vector(); + m.merged = std::vector(); } if (input.segments.size() < 2) logger.Warning("Only {} sweep goes into the joint merge", input.segments.size());