diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 2d33580de..abd557966 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -462,7 +462,8 @@ in the presence of an indexing ambiguity" (2018), Acta Cryst. D74, 405-410 [doi:10.1107/S2059798318002978](https://doi.org/10.1107/S2059798318002978). **Reference sweep of a joint merge** — the sweep the others are put on the axes of is the one whose own -merge reaches the finest resolution, as KAMO takes the data set of the smallest Wilson B as its reference: +merge reaches the finest resolution (among the sweeps with a measured ISa), as KAMO takes the data set of +the smallest Wilson B as its reference: K. Yamashita, K. Hirata and M. Yamamoto, "KAMO: towards automated data processing for microcrystals" (2018), Acta Cryst. D74, 441-449 [doi:10.1107/S2059798318004576](https://doi.org/10.1107/S2059798318004576). diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index afbb9fbe6..9c8724485 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -373,7 +373,7 @@ An internal-consistency term can inflate uncertainties when multiple observation **Outlier rejection.** Before merging, an observation farther than $N\sigma$ from its reflection's median is rejected (`--reject-outliers`; default 6 on rotation data and 10 on stills, whose partials scatter with heavier tails; 0 turns it off). -**Several sweeps.** Rotation sweeps given together are one dataset by default (`--sweeps joint`): each is scaled on its own, the sweeps are put on the axes of a reference sweep (the one reaching the finest resolution, ties broken by the higher ISa) and merged in two levels, as XSCALE merges separately scaled inputs; `--sweeps separate` processes each alone. See [RUGNUX_CLI](RUGNUX_CLI.md#inputs). +**Several sweeps.** Rotation sweeps given together are one dataset by default (`--sweeps joint`): each is scaled on its own, the sweeps are put on the axes of a reference sweep (the one reaching the finest resolution among those with a measured ISa, ties broken by the higher ISa) and merged in two levels, as XSCALE merges separately scaled inputs; `--sweeps separate` processes each alone. See [RUGNUX_CLI](RUGNUX_CLI.md#inputs). ### 10.5 Merging statistics diff --git a/docs/RUGNUX_CLI.md b/docs/RUGNUX_CLI.md index 6718cd50d..b3ec02335 100644 --- a/docs/RUGNUX_CLI.md +++ b/docs/RUGNUX_CLI.md @@ -107,14 +107,14 @@ so are sweeps at different wavelengths, `-s`/`-e` and `--force-still` with sever By default (`--sweeps joint`) the sweeps are one dataset, merged in two levels as XSCALE merges the data sets XDS made of each sweep: each sweep is processed on its own first - indexed, integrated, its partials combined and scaled, its space group determined - writing only its log `_s_log.txt`; -the sweeps are then put on the axes of the best of them (the finest resolution of its own merge, ties: the higher ISa) - the change +the sweeps are then put on the axes of the best of them (the finest resolution of its own merge among the sweeps whose +error model resolved an ISa, ties: the higher ISa) - the change of basis, an axis permutation or a twin law of the lattice, under which their merged intensities correlate best; where the lattice allows more than one indexing, which one each sweep takes is settled over all pairs of sweeps at once - and each is judged by its isomorphism ratio with it (see [RUGNUX_REPORT](RUGNUX_REPORT.md)): clearly about 1 it joins, clearly below it (another crystal) it is left out, and in between the merge decides - it -joins unless it hurts - adding it raises the held-out residual of the others' merge by at most 1% (two -standard errors up), or it adds at least 1% to their unique reflections while that residual is not worse -beyond its noise (a sweep left out is named in the log). The crystal has one cell: a joined sweep whose own cell is more than 0.3% off that of +joins unless it is shown to hurt: adding it raises the held-out residual of the others' merge by more than 1%, +two standard errors down (a sweep left out is named in the log). The crystal has one cell: a joined sweep whose own cell is more than 0.3% off that of the joined sweep reaching the finest resolution (where the cell is measured rather than the distance) is processed a second time, from the distance that would make the two cells agree, and that run replaces the first where its intensities agree better with that sweep's (a lower held-out residual). Their fulls are diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index ac86760bc..47d07a21f 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -331,7 +331,8 @@ A run of several sweeps merged as one (`--sweeps joint`, the default with severa **1b. SWEEPS**; a run of one input has no such section and no such keys. - **`SWEEP_COUNT=`** the sweeps given, **`SWEEP_REFERENCE=`** the one (1-based) the others were put on - the axes of - the one reaching the finest resolution, ties broken by the higher ISa - and **`SWEEP_DECAY_MODEL=`** how the decay over the + the axes of - the one reaching the finest resolution among the sweeps whose error model resolved an ISa, ties + broken by the higher ISa - and **`SWEEP_DECAY_MODEL=`** how the decay over the sweeps was modelled: `dose-linear` (a B rising along the dose order), `per-sweep B` (one B a sweep) or `none`, whichever predicted held-out frames better. - Per sweep `k`: **`SWEEP_k_INPUT=`**, **`SWEEP_k_IMAGES=`**, **`SWEEP_k_INCLUDED=`**; what the sweep @@ -342,7 +343,7 @@ A run of several sweeps merged as one (`--sweeps joint`, the default with severa own CC1/2 (CC_AB / (CC*_A CC*_B), Karplus & Diederichs 2012), on the reflections both have in the coarser half of their common range, with a bootstrap error - 0.997-1.001 for sweeps of one crystal or isomorphous crystals, 0.95-0.96 for non-isomorphous ones. **`SWEEP_k_DECISION=`** says what became of the sweep: - `reference` (the sweep whose own merge reaches the finest resolution, ties: the higher ISa, whose axes and cell the joint merge takes), + `reference` (the sweep whose own merge reaches the finest resolution among those with an ISa, ties: the higher ISa, whose axes and cell the joint merge takes), `isomorphous` (ratio - 2 se >= 0.99: joined), `not isomorphous` (ratio + 2 se < 0.97: left out), `doubtful: kept` / `doubtful: left out` (in between, or too few common reflections for a ratio: decided on the merge), `no axes` (it could not be put on the reference's axes) or `no merge`. For a doubtful sweep, @@ -350,8 +351,7 @@ A run of several sweeps merged as one (`--sweeps joint`, the default with severa residual of the merge of the joined sweeps - each half of their observations predicted from the other half, without and with this sweep's added, over the reflections it has, after a scale and B factor onto the reference - with its bootstrap error, and **`SWEEP_k_HELD_OUT_ADDED=`**, the unique reflections it adds to - that merge; the sweep joins where the change is at most 1% two standard errors up, or it adds at least 1% - to the merge's unique reflections and the change is not above twice its error; **`SWEEP_k_DISTANCE=`**, the detector + that merge; the sweep joins unless the change is above 1% two standard errors down; **`SWEEP_k_DISTANCE=`**, the detector distance its run ended at, and **`SWEEP_k_REPROCESSED_FROM=`** where the sweep was processed a second time (see [RUGNUX_CLI](RUGNUX_CLI.md)): the distance its first run ended at; and, for a sweep in the joint merge, **`SWEEP_k_REINDEX_OP=`** (onto the reference's axes) with diff --git a/rugnux/RugnuxJointSweeps.cpp b/rugnux/RugnuxJointSweeps.cpp index 14f45f200..6da5805fd 100644 --- a/rugnux/RugnuxJointSweeps.cpp +++ b/rugnux/RugnuxJointSweeps.cpp @@ -55,6 +55,7 @@ SweepIsomorphism MeasureSweepIsomorphism(const std::vector &a, const sa.w[0] > 0 && sa.w[1] > 0, sb.w[0] > 0 && sb.w[1] > 0, sa.d, k}); } SweepIsomorphism out; + out.n = static_cast(rows.size()); constexpr size_t MIN_ROWS = 200; // 100 in the half the ratio is measured on if (rows.size() < MIN_ROWS) return out; @@ -63,6 +64,7 @@ SweepIsomorphism MeasureSweepIsomorphism(const std::vector &a, const return x.d != y.d ? x.d > y.d : x.key < y.key; // ties by index, not by the map's hash order }); rows.resize(rows.size() / 2); + out.n = static_cast(rows.size()); const auto ratio = [&](const std::vector &idx) { const auto cc = [&](auto x, auto y, auto use) { @@ -88,7 +90,6 @@ SweepIsomorphism MeasureSweepIsomorphism(const std::vector &a, const std::vector idx(rows.size()); std::iota(idx.begin(), idx.end(), 0); out.ratio = ratio(idx); - out.n = static_cast(rows.size()); constexpr int BOOTSTRAP = 200; std::mt19937 rng(12345); std::uniform_int_distribution pick(0, rows.size() - 1); @@ -366,9 +367,12 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps out.sweeps.resize(sweeps.size()); if (sweeps.empty()) return out; - // The best sweep on its own: the one reaching the finest resolution, ties broken by the higher ISa (a - // sweep that merged nothing comes last). + // The best sweep on its own: of the sweeps whose error model resolved an ISa, the one reaching the + // finest resolution, ties broken by the higher ISa (a sweep that merged nothing comes last). + const auto accuracy_known = [&](size_t k) { return sweeps[k].isa > 0.0; }; // false on NaN const auto better = [&](size_t k, size_t than) { + if (accuracy_known(k) != accuracy_known(than)) + return accuracy_known(k); const double dk = std::isfinite(sweeps[k].d_min) ? sweeps[k].d_min : INFINITY; const double dt = std::isfinite(sweeps[than].d_min) ? sweeps[than].d_min : INFINITY; if (dk != dt) @@ -413,6 +417,12 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps JointSweepRelation::Sweep &s = out.sweeps[k]; if (static_cast(k) == out.reference || !s.decision.empty()) continue; + if (s.cc < MIN_REFERENCE_CC) { + s.decision = "no axes"; + s.why_left_out = fmt::format("its intensities correlate with the reference sweep's on no axes " + "(best CC {:.3f})", s.cc); + continue; + } const SweepIsomorphism iso = MeasureSweepIsomorphism(sweeps[k].fulls, s.cob, ref.fulls, gemmi::Op::identity(), ref.space_group); s.iso_ratio = iso.ratio; @@ -433,11 +443,6 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps iso.ratio, iso.se, iso.n, LEAVE_RATIO); continue; } - } else if (s.cc < MIN_REFERENCE_CC) { - s.decision = "no axes"; - s.why_left_out = fmt::format("its intensities correlate with the reference sweep's on no axes " - "(best CC {:.3f})", s.cc); - continue; } doubtful.push_back(k); } @@ -457,21 +462,14 @@ JointSweepRelation RelateJointSweeps(const std::vector &sweeps s.held_out_se = g.se; s.held_out_n = g.n; s.held_out_added = g.added; - s.included = std::isfinite(g.change) - && (g.change + 2.0 * g.se <= MIN_HELD_OUT_DEGRADATION - || (static_cast(g.added) >= MIN_ADDED_FRACTION * static_cast(g.base_unique) - && g.added > 0 && std::isfinite(g.se) && g.change <= 2.0 * g.se)); + s.included = !(std::isfinite(g.change) && std::isfinite(g.se) + && g.change - 2.0 * g.se > MIN_HELD_OUT_DEGRADATION); s.decision = s.included ? "doubtful: kept" : "doubtful: left out"; - if (s.included) - continue; - if (std::isfinite(g.change)) + if (!s.included) s.why_left_out = fmt::format("its isomorphism ratio ({:.3f} +- {:.3f}) leaves it in doubt, and adding " "it raises the held-out residual of the merge ({:+.1f}% +- {:.1f}% on " "{} reflections, {} unique reflections added)", s.iso_ratio, s.iso_se, 100.0 * g.change, 100.0 * g.se, g.n, g.added); - else - s.why_left_out = fmt::format("its isomorphism ratio ({:.3f}) leaves it in doubt, and it shares too " - "few reflections with the merge to test it there", s.iso_ratio); } return out; } diff --git a/rugnux/RugnuxJointSweeps.h b/rugnux/RugnuxJointSweeps.h index 0acc26ab7..4150100c9 100644 --- a/rugnux/RugnuxJointSweeps.h +++ b/rugnux/RugnuxJointSweeps.h @@ -76,40 +76,43 @@ HeldOutGain MeasureHeldOutGain(const std::vector * const std::vector &d, const gemmi::Op &cob_d, 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 sweep -// whose own merge reaches the finest resolution (ties: the higher ISa) - the best-diffracting, least -// damaged one, as KAMO (Yamashita et al. (2018) Acta Cryst. D74, 441-449) picks the data set of the smallest Wilson B - not the largest: a long but poor -// sweep is a poor yardstick for the others, and one that stays at low resolution shares too few -// reflections with the rest. Each other sweep is put on its axes (the change of basis, an -// axis permutation or a twin law of the lattice, ChooseReferenceSetting, under which its merged +// 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 +// resolved an ISa, the one whose own merge reaches the finest resolution (ties: the higher ISa) - the +// least damaged one, as KAMO (Yamashita et al. (2018) Acta Cryst. D74, 441-449) picks the data set of the +// smallest Wilson B - not the largest: a long but poor sweep is a poor yardstick for the others, and one +// that stays at low resolution shares too few reflections with the rest. The resolution alone is not +// enough: a lattice that integrated a few hundred reflections of noise out to the edge of the detector +// reaches the finest resolution of all, with 2% of the data merged and no ISa its error model could +// resolve, and the merge made on its axes held nothing. Nor is a sweep's own CC1/2 or its unique count +// a better merit (both were tried): CC1/2 saturates on a small cell and reads precision, which a long +// sweep of a poor crystal has, and a count rewards a sweep merged in too high a symmetry; both moved the +// reference off the sweep the others agree with. Each other sweep is put on the reference's axes (the +// change of basis, an axis permutation or a twin law of the lattice, ChooseReferenceSetting, under which its merged // intensities correlate best with the reference's) and judged by its isomorphism ratio with it // (SweepIsomorphism), by effect size, not significance: two sweeps of one crystal differ systematically by // 0.001-0.003 (absorption, dose), far beyond the ratio's bootstrap error over 10^4 reflections, so a test // against 1 or against the reference's own halves rejects true pairs. Ratio - 2 se at or above JOIN_RATIO: // the same intensities, it joins; ratio + 2 se below LEAVE_RATIO: another crystal, it is left out. In // between, or where too few reflections are common to measure the ratio, the sweep is doubtful and the -// merge decides: it joins where adding it raises the held-out residual of the clearly joined sweeps -// (MeasureHeldOutGain) by no more than MIN_HELD_OUT_DEGRADATION, two standard errors up - more data of the -// same crystal is worth having unless it costs something that matters - or where it adds unique -// reflections to them while that residual is not worse beyond twice its error (a weak or damaged but -// isomorphous wedge still adds coverage the residual cannot see). A sweep that cannot be put on the -// reference's axes (too few common reflections, or a best correlation below MIN_REFERENCE_CC) is left out. +// merge decides: it joins unless adding it is shown to raise the held-out residual of the clearly joined +// sweeps (MeasureHeldOutGain) by more than MIN_HELD_OUT_DEGRADATION, two standard errors down. A doubtful +// sweep already correlates with the reference on its axes and the ratio could not tell it from the same +// crystal, so more of its data is worth having unless the merge measures a cost; one whose effect the +// merge cannot measure - a cell of a few dozen unique reflections, where no statistic over reflections +// resolves anything, or too few reflections in common for the residual at all - joins on that +// correlation. A sweep that cannot be put on the reference's axes (too few common reflections, or a best +// correlation below MIN_REFERENCE_CC) is left out. // // The band edges come from the pairs measured when this was written: sweeps of one crystal 0.998-1.001 // (three pairs from four sweeps of a cubic crystal, one tetragonal pair), two isomorphous cubic crystals // 0.997 +- 0.002; two non-isomorphous tetragonal crystals 0.951 +- 0.003, two non-isomorphous trigonal // crystals 0.962 +- 0.002. constexpr double JOIN_RATIO = 0.99; -// What "adds unique reflections" means for a doubtful sweep: at least 1% more than the joined sweeps have. -// Below that its coverage changes the completeness by less than the last digit a report prints, so it -// buys nothing a user would see, and a handful of reflections at the edge of its range must not carry a -// doubtful sweep in. The small wedges this is for add 10-60% (a 50 deg wedge to two others: ~30%). -constexpr double MIN_ADDED_FRACTION = 0.01; -// The rise of the held-out residual, two standard errors up, a doubtful sweep may cause and still join: one -// value for every run. Measured, sweeps of the same crystal raised it by 0.1-0.2% +- 0.03-0.05% - beyond -// twice its error or not depending on a small change of integration that left their own merges as good or -// better - and a change that size is no reason to lose their data. Two standard errors up, so that a sweep -// whose effect is not measured at all (a few reflections in common) does not pass for a harmless one. +// The rise of the held-out residual a doubtful sweep may cause and still join: one value for every run. +// Measured, sweeps of the same crystal raised it by 0.1-0.2% +- 0.03-0.05% - beyond twice its error or +// not depending on a small change of integration that left their own merges as good or better - and a +// change that size is no reason to lose their data. constexpr double MIN_HELD_OUT_DEGRADATION = 0.01; constexpr double LEAVE_RATIO = 0.97; diff --git a/tests/RugnuxJointSweepsTest.cpp b/tests/RugnuxJointSweepsTest.cpp index 2d67542be..6fd660ee7 100644 --- a/tests/RugnuxJointSweepsTest.cpp +++ b/tests/RugnuxJointSweepsTest.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -33,6 +34,7 @@ JointMemberSweep MakeSweep(const std::map, double> &truth, co s.n_frames = n_frames; s.cell = UnitCell{A, A, C, 90, 90, 90}; s.space_group = P4(); + s.isa = 1000.0 / sigma; // what its error model would resolve: the noise level against the intensities for (const auto &[hkl, I]: truth) { const auto m = indexed_as.apply_to_hkl({{hkl[0], hkl[1], hkl[2]}}); const auto d = static_cast(1.0 / std::sqrt((hkl[0] * hkl[0] + hkl[1] * hkl[1]) / (A * A) @@ -55,12 +57,26 @@ JointMemberSweep MakeSweep(const std::map, double> &truth, co r.h = m[0]; r.k = m[1]; r.l = m[2]; r.I = static_cast(sum / 4); r.sigma = static_cast(sigma / 2); + r.I_half[0] = (s.fulls[s.fulls.size() - 4].I + s.fulls[s.fulls.size() - 2].I) / 2; + r.I_half[1] = (s.fulls[s.fulls.size() - 3].I + s.fulls[s.fulls.size() - 1].I) / 2; r.d = d; s.merged.push_back(r); } + s.d_min = s.merged.empty() ? NAN : std::min_element(s.merged.begin(), s.merged.end(), [](const auto &a, const auto &b) { + return a.d < b.d; + })->d; return s; } +// The reflections of `truth` finer than `d` dropped: a sweep that stops there. +std::map, double> CoarserThan(const std::map, double> &truth, double d) { + std::map, double> t; + for (const auto &[hkl, I]: truth) + if (1.0 / std::sqrt((hkl[0] * hkl[0] + hkl[1] * hkl[1]) / (A * A) + hkl[2] * hkl[2] / (C * C)) >= d) + t[hkl] = I; + return t; +} + std::map, double> Truth(uint32_t seed) { std::mt19937 rng(seed); std::exponential_distribution wilson(1.0); @@ -81,12 +97,9 @@ TEST_CASE("RelateJointSweeps_AmbiguityRecoveredAndAStrangerLeftOut", "[Rugnux]") const auto truth = Truth(1); const gemmi::Op twin = gemmi::parse_triplet("y,x,-z"); // h,k,l -> k,h,-l std::vector sweeps; - sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 2)); - sweeps.push_back(MakeSweep(truth, twin, 80, 3)); - sweeps.push_back(MakeSweep(Truth(99), gemmi::Op::identity(), 50, 4)); // another crystal - sweeps[0].isa = 20.0; // the reference: the best ISa - sweeps[1].isa = 10.0; - sweeps[2].isa = 15.0; + sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 2, 10.0)); // the reference: the best measured + sweeps.push_back(MakeSweep(truth, twin, 80, 3, 20.0)); + sweeps.push_back(MakeSweep(Truth(99), gemmi::Op::identity(), 50, 4, 20.0)); // another crystal const auto relation = RelateJointSweeps(sweeps); REQUIRE(relation.reference == 0); @@ -98,8 +111,9 @@ TEST_CASE("RelateJointSweeps_AmbiguityRecoveredAndAStrangerLeftOut", "[Rugnux]") const auto back = relation.sweeps[1].cob.apply_to_hkl(twin.apply_to_hkl({{3, 1, 2}})); CHECK(P4Key(back[0], back[1], back[2]) == P4Key(3, 1, 2)); CHECK(relation.sweeps[1].decision == "isomorphous"); + // another crystal's intensities correlate with the reference's on no axes CHECK_FALSE(relation.sweeps[2].included); - CHECK(relation.sweeps[2].decision == "not isomorphous"); + CHECK(relation.sweeps[2].decision == "no axes"); const auto input = BuildJointMergeInput(sweeps, relation); REQUIRE(input.segments.size() == 2); @@ -120,8 +134,8 @@ TEST_CASE("RelateJointSweeps_AmbiguityRecoveredAndAStrangerLeftOut", "[Rugnux]") } TEST_CASE("RelateJointSweeps_DoubtfulSweepDecidedOnTheMerge", "[Rugnux]") { - // Sweep 1: the same crystal measured noisily, its ratio uncertain; sweep 2: precise data of a crystal - // whose intensities differ by a few percent from the reference's. Both fall in the doubtful band, and + // Sweep 1: the same crystal measured noisily, its ratio uncertain; sweep 2: a crystal whose + // intensities differ by a few percent from the reference's, measured well enough to show it. Both fall in the doubtful band, and // the held-out residual keeps the first and leaves out the second. const auto truth = Truth(1); auto shifted = truth; @@ -130,10 +144,9 @@ TEST_CASE("RelateJointSweeps_DoubtfulSweepDecidedOnTheMerge", "[Rugnux]") { for (auto &[hkl, I]: shifted) I *= std::exp(change(rng)); std::vector sweeps; - sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 2, 100.0)); + sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 2, 50.0)); sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 3, 1500.0)); - sweeps.push_back(MakeSweep(shifted, gemmi::Op::identity(), 100, 4, 10.0)); - sweeps[0].isa = 20.0; + sweeps.push_back(MakeSweep(shifted, gemmi::Op::identity(), 100, 4, 100.0)); const auto relation = RelateJointSweeps(sweeps); REQUIRE(relation.reference == 0); @@ -146,3 +159,49 @@ TEST_CASE("RelateJointSweeps_DoubtfulSweepDecidedOnTheMerge", "[Rugnux]") { CHECK_FALSE(relation.sweeps[2].included); CHECK(relation.sweeps[2].held_out_change > 0.0); } + +TEST_CASE("RelateJointSweeps_ReferenceIsTheBestMeasuredSweepNotTheFinestResolution", "[Rugnux]") { + // Sweep 1: the crystal measured well to 5 A. Sweep 2: the same reflections, and some finer ones, all of + // them noise - what a lattice that integrated nothing but background to the edge of the detector merges, + // with no ISa its error model could resolve. Its own merge reaches the finest resolution; the sweep + // with a measured accuracy is the reference. + const auto truth = Truth(1); + std::vector sweeps; + sweeps.push_back(MakeSweep(CoarserThan(truth, 5.0), gemmi::Op::identity(), 100, 2, 10.0)); + sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 3, 5000.0)); + sweeps[1].isa = NAN; + REQUIRE(sweeps[1].d_min < sweeps[0].d_min); + + const auto relation = RelateJointSweeps(sweeps); + INFO("sweep 2: cc " << relation.sweeps[1].cc << ", ratio " << relation.sweeps[1].iso_ratio << " +- " + << relation.sweeps[1].iso_se << ", held-out " << relation.sweeps[1].held_out_change << " +- " + << relation.sweeps[1].held_out_se << ", " << relation.sweeps[1].decision); + CHECK(relation.reference == 0); + CHECK(relation.sweeps[0].decision == "reference"); + CHECK_FALSE(relation.sweeps[1].included); + CHECK(relation.sweeps[1].decision == "no axes"); +} + +TEST_CASE("RelateJointSweeps_ASweepTheMergeCannotMeasureJoins", "[Rugnux]") { + // A cell of a few dozen unique reflections: too few for the isomorphism ratio, and the held-out + // residual of the merge moves by less than its own noise. Nothing measured refuses the sweep, so it + // joins - the second measurement of the same reflections is what a small-molecule merge needs. + std::map, double> small; + int n = 0; + for (const auto &[hkl, I]: Truth(1)) + if (n++ % 50 == 0) + small[hkl] = I; + REQUIRE(small.size() < 200); + std::vector sweeps; + sweeps.push_back(MakeSweep(small, gemmi::Op::identity(), 100, 2, 200.0)); + sweeps.push_back(MakeSweep(small, gemmi::Op::identity(), 100, 3, 300.0)); + + const auto relation = RelateJointSweeps(sweeps); + REQUIRE(relation.reference == 0); + INFO("held-out change " << relation.sweeps[1].held_out_change << " +- " << relation.sweeps[1].held_out_se + << " on " << relation.sweeps[1].held_out_n << " reflections"); + CHECK(std::isnan(relation.sweeps[1].iso_ratio)); + CHECK(std::isnan(relation.sweeps[1].held_out_change)); + CHECK(relation.sweeps[1].decision == "doubtful: kept"); + CHECK(relation.sweeps[1].included); +}