Merge branch 'rc186-fix-multisweep' into rc186-merge2

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<MergedReflectio
                   [&](int a, int b) { return merged[a].d != merged[b].d ? merged[a].d > merged[b].d : a < b; });
         const int n_beyond = static_cast<int>(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<int>(kMinModulationBand, static_cast<int>(mod_idx.size()))));
+        mod_idx.resize(std::min(kMaxModulationBand, std::max(n_beyond, std::min<int>(kMinModulationBand,
+                                                                                     static_cast<int>(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<MergedReflectio
                     cls_hkl.push_back({merged[i].h, merged[i].k, merged[i].l});
                     cls_e2.push_back(e2_of[i]);
                 }
+            // The ladder up to the lowest kMaxClassReflections orders.
+            if (cls_hkl.size() > 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

rugnux/RugnuxJointSweeps.cpp

@@ -5,12 +5,14 @@
 
 #include <algorithm>
 #include <cmath>
+#include <future>
 #include <map>
 #include <numeric>
 #include <random>
 #include <unordered_map>
 
 #include "../common/CrystalLattice.h"
+#include "../common/Logger.h"
 #include <spdlog/fmt/fmt.h>
 
 #include "../image_analysis/scale_merge/HKLKey.h"
@@ -109,10 +111,11 @@ SweepIsomorphism MeasureSweepIsomorphism(const std::vector<ScaledFull> &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<const std::vector<ScaledFull> *> &base,
-                               const std::vector<gemmi::Op> &base_cob,
-                               const std::vector<ScaledFull> &d, const gemmi::Op &cob_d,
-                               const gemmi::SpaceGroup &space_group) {
+std::vector<HeldOutGain> MeasureHeldOutGains(const std::vector<const std::vector<ScaledFull> *> &base,
+                                              const std::vector<gemmi::Op> &base_cob,
+                                              const std::vector<const std::vector<ScaledFull> *> &ds,
+                                              const std::vector<gemmi::Op> &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<ScaledFull> &fulls, const gemmi::Op &cob) {
@@ -180,84 +183,102 @@ HeldOutGain MeasureHeldOutGain(const std::vector<const std::vector<ScaledFull> *
         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<ScaledFull> &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<Row> 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<double>(o.sigma) * o.sigma);
-            wid += o.I / (static_cast<double>(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<double>(obs[i].sigma) * obs[i].sigma);
-                wi += obs[i].I / (static_cast<double>(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<Row> 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<double>(o.sigma) * o.sigma);
+                wid += o.I / (static_cast<double>(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<double>(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<double>(obs[i].sigma) * obs[i].sigma);
+                    wi += obs[i].I / (static_cast<double>(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<double>(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<int64_t>(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<size_t> &idx) {
-        std::vector<double> 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<double> &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<int64_t>(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<size_t> &idx) {
+            std::vector<double> 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<double> &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<size_t> idx(rows.size());
+        std::iota(idx.begin(), idx.end(), 0);
+        out.change = change(idx);
+        out.n = static_cast<int64_t>(rows.size());
+        constexpr int BOOTSTRAP = 200;
+        std::mt19937 rng(12345);
+        std::uniform_int_distribution<size_t> 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<size_t> idx(rows.size());
-    std::iota(idx.begin(), idx.end(), 0);
-    out.change = change(idx);
-    out.n = static_cast<int64_t>(rows.size());
-    constexpr int BOOTSTRAP = 200;
-    std::mt19937 rng(12345);
-    std::uniform_int_distribution<size_t> 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<std::future<HeldOutGain>> 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<HeldOutGain> out;
+    for (auto &g : gains)
+        out.push_back(g.get());
     return out;
 }
 
+HeldOutGain MeasureHeldOutGain(const std::vector<const std::vector<ScaledFull> *> &base,
+                               const std::vector<gemmi::Op> &base_cob,
+                               const std::vector<ScaledFull> &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<JointMemberSweep> &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<int>(k) == out.reference)
             continue;
@@ -410,9 +433,14 @@ JointSweepRelation RelateJointSweeps(const std::vector<JointMemberSweep> &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<size_t> doubtful;
+    TimingMark iso_mark(logger, "joint sweeps: isomorphism");
+    // Each sweep's ratio is measured on its own, all of them at once.
+    std::vector<std::future<SweepIsomorphism>> isomorphism(sweeps.size());
     for (size_t k = 0; k < sweeps.size(); ++k) {
         JointSweepRelation::Sweep &s = out.sweeps[k];
         if (static_cast<int>(k) == out.reference || !s.decision.empty())
@@ -423,8 +451,16 @@ JointSweepRelation RelateJointSweeps(const std::vector<JointMemberSweep> &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<size_t> 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<JointMemberSweep> &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<const std::vector<ScaledFull> *> base{&ref.fulls};
     std::vector<gemmi::Op> base_cob{gemmi::Op::identity()};
     for (size_t k = 0; k < sweeps.size(); ++k)
@@ -455,9 +494,18 @@ JointSweepRelation RelateJointSweeps(const std::vector<JointMemberSweep> &sweeps
             base.push_back(&sweeps[k].fulls);
             base_cob.push_back(out.sweeps[k].cob);
         }
+    std::vector<const std::vector<ScaledFull> *> doubtful_fulls;
+    std::vector<gemmi::Op> 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<HeldOutGain> 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<double>(CrystalLattice(cell).CalcVolume()));
 }
 
-HeldOutGain HeldOutAgainst(const std::vector<const JointMemberSweep *> &base, const std::vector<gemmi::Op> &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<HeldOutGain> HeldOutAgainst(const std::vector<const JointMemberSweep *> &base,
+                                        const std::vector<gemmi::Op> &base_cob, const JointMemberSweep &reference,
+                                        const std::vector<const JointMemberSweep *> &ds) {
     std::vector<const std::vector<ScaledFull> *> 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<ScaledFull> none;
+    std::vector<const std::vector<ScaledFull> *> d_fulls;
+    std::vector<gemmi::Op> 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<JointMemberSweep> &sweeps, const JointSweepRelation &relation) {
@@ -518,11 +574,24 @@ JointMergeInput BuildJointMergeInput(const std::vector<JointMemberSweep> &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<int64_t>(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<size_t> 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<int64_t>(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

rugnux/RugnuxJointSweeps.h

@@ -75,6 +75,12 @@ HeldOutGain MeasureHeldOutGain(const std::vector<const std::vector<ScaledFull> *
                                const std::vector<gemmi::Op> &base_cob,
                                const std::vector<ScaledFull> &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<HeldOutGain> MeasureHeldOutGains(const std::vector<const std::vector<ScaledFull> *> &base,
+                                             const std::vector<gemmi::Op> &base_cob,
+                                             const std::vector<const std::vector<ScaledFull> *> &ds,
+                                             const std::vector<gemmi::Op> &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<JointMemberSweep> &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<const JointMemberSweep *> &base, const std::vector<gemmi::Op> &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<HeldOutGain> HeldOutAgainst(const std::vector<const JointMemberSweep *> &base,
+                                        const std::vector<gemmi::Op> &base_cob, const JointMemberSweep &reference,
+                                        const std::vector<const JointMemberSweep *> &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

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();

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<float> 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<float> distance_mm, std::optional<UnitCell> 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<std::string> 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<float>(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<HeldOutGain> 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<ScaledFull>();
+        m.merged = std::vector<MergedReflection>();
     }
     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;