Joint sweeps: no reference without a measured ISa, and a doubtful sweep joins unless shown to hurt
The reference sweep was the one whose own merge reached the finest resolution. A lattice that integrated a few hundred reflections of noise out to the detector edge reaches the finest resolution of all with 2% of the data merged, became the reference, and the real sweeps were then judged against it: they shared too few reflections with it and were left out, and the joint merge held one lattice of noise. Such a lattice's error model resolves no ISa, and that is the whole change to the rule: of the sweeps with a measured ISa, the one reaching the finest resolution, ties by ISa. Two other merits were tried against the multi-sweep rows and rejected: a sweep's own CC1/2 saturates on a small cell (the smallest of seven kappa runs won by noise, and nothing could be measured against it) and reads precision (a long sweep of a poor crystal won over the best one); coverage times ISa rewards a sweep merged in too high a symmetry and a minority twin domain with the higher ISa. A doubtful sweep (isomorphism ratio in the band, or not measurable) was kept only where the held-out residual of the merge was shown not to rise (two standard errors up), or where it added coverage. On a cell of a few dozen unique reflections no statistic over reflections resolves anything, so a second sweep of the same crystal was left out on a residual change of +1.6% +- 1.3% - and a merge of both was clearly the better one. It now joins unless adding it is shown to raise the residual by more than the allowance, two standard errors down; one the merge cannot measure at all joins on its correlation with the reference. The added-coverage clause is subsumed and goes. On seven kappa runs of a small molecule this joins one run more than before, with a better refinement R1. A sweep whose intensities correlate with the reference's on no axes (best CC below MIN_REFERENCE_CC) is left out before the isomorphism test, as the header always said; it used to reach the band on a ratio measured by chance. The ratio's reflection count is now reported where there were too few to measure it, instead of 0. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
This commit is contained in:
1 parent
22dd8d8528
commit
0182384f42
7 files changed
+123
-62
No files matched your search
@@ -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).
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
+4
-4
@@ -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 `<prefix>_s<k>_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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -55,6 +55,7 @@ SweepIsomorphism MeasureSweepIsomorphism(const std::vector<ScaledFull> &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<int64_t>(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<ScaledFull> &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<int64_t>(rows.size());
|
||||
|
||||
const auto ratio = [&](const std::vector<size_t> &idx) {
|
||||
const auto cc = [&](auto x, auto y, auto use) {
|
||||
@@ -88,7 +90,6 @@ SweepIsomorphism MeasureSweepIsomorphism(const std::vector<ScaledFull> &a, const
|
||||
std::vector<size_t> idx(rows.size());
|
||||
std::iota(idx.begin(), idx.end(), 0);
|
||||
out.ratio = ratio(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);
|
||||
@@ -366,9 +367,12 @@ JointSweepRelation RelateJointSweeps(const std::vector<JointMemberSweep> &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<JointMemberSweep> &sweeps
|
||||
JointSweepRelation::Sweep &s = out.sweeps[k];
|
||||
if (static_cast<int>(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<JointMemberSweep> &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<JointMemberSweep> &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<double>(g.added) >= MIN_ADDED_FRACTION * static_cast<double>(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;
|
||||
}
|
||||
|
||||
+25
-22
@@ -76,40 +76,43 @@ HeldOutGain MeasureHeldOutGain(const std::vector<const std::vector<ScaledFull> *
|
||||
const std::vector<ScaledFull> &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;
|
||||
|
||||
|
||||
@@ -4,6 +4,7 @@
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
#include <map>
|
||||
#include <random>
|
||||
|
||||
@@ -33,6 +34,7 @@ JointMemberSweep MakeSweep(const std::map<std::array<int, 3>, 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<float>(1.0 / std::sqrt((hkl[0] * hkl[0] + hkl[1] * hkl[1]) / (A * A)
|
||||
@@ -55,12 +57,26 @@ JointMemberSweep MakeSweep(const std::map<std::array<int, 3>, double> &truth, co
|
||||
r.h = m[0]; r.k = m[1]; r.l = m[2];
|
||||
r.I = static_cast<float>(sum / 4);
|
||||
r.sigma = static_cast<float>(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<std::array<int, 3>, double> CoarserThan(const std::map<std::array<int, 3>, double> &truth, double d) {
|
||||
std::map<std::array<int, 3>, 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<std::array<int, 3>, double> Truth(uint32_t seed) {
|
||||
std::mt19937 rng(seed);
|
||||
std::exponential_distribution<double> 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<JointMemberSweep> 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<JointMemberSweep> 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<JointMemberSweep> 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<std::array<int, 3>, 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<JointMemberSweep> 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);
|
||||
}
|
||||
Reference in new issue
Block a user