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/RugnuxJointSweeps.cpp b/rugnux/RugnuxJointSweeps.cpp index 6da5805fd..db0a59943 100644 --- a/rugnux/RugnuxJointSweeps.cpp +++ b/rugnux/RugnuxJointSweeps.cpp @@ -5,12 +5,14 @@ #include #include +#include #include #include #include #include #include "../common/CrystalLattice.h" +#include "../common/Logger.h" #include #include "../image_analysis/scale_merge/HKLKey.h" @@ -109,10 +111,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 +183,102 @@ 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))); + // ... 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) + 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; } +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 +411,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,9 +433,14 @@ 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(); - std::vector doubtful; + TimingMark iso_mark(logger, "joint sweeps: isomorphism"); + // 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()) @@ -423,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; @@ -447,7 +483,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 +494,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 +526,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) { @@ -518,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/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/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 98ad0d0c7..b366aa678 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 fe6dea69c..9ca68f76a 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(); @@ -4094,13 +4103,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; @@ -4151,16 +4163,36 @@ int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { std::vector inputs; 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, "")); + members.push_back(run_member(k, std::nullopt, std::nullopt, "")); inputs.push_back(plan.sweeps[k].input); } // 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 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 (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 +4206,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 +4218,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 +4229,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; @@ -4228,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()); @@ -4248,6 +4283,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;