leonarski_fandClaude Opus 5.5
439ca10db1
A still indexed on few spots often gets a lattice that explains them as well as the true one: an axis multiple, whose reciprocal lattice holds every true node plus the nodes between, or the true lattice in another setting. Its cell is then not the run's, the merge leaves it out, and its reflections at the true nodes - as good as any - are lost. On two 25k-frame XFEL runs of a tetragonal crystal these were 26-30% of the indexed frames; for those holding the true cell as a sublattice, the reflections at true nodes had the I/sigma of accepted frames and those between them none. RelabelOnSublattice (new, scale_merge/SublatticeRelabel) compares lattices, not cell parameters: it looks for integer combinations of the frame's reciprocal axes whose full metric is the run's cell within the merge's own tolerances, volume ratio 1-4 (another setting, or a 2x-4x supercell). Where several such bases exist - a nearly symmetric supercell fits along more than one axis - the one with the strongest mean I/sigma at its nodes is taken, since only the true lattice has signal there. The frame's indices are re-derived on that basis (h = h' M^-1), reflections that are not nodes are dropped, as are those finer than any frame already on the run's cell reaches (the wrong lattice was also predicted to a wrong resolution), and the frame gets the run's basis in its own orientation. It runs before the first stills merge, and in --mode scale before the resolutions are recomputed from the cell. --mode scale, both runs (indicative; the mx pairs are queued): 1060 and 920 crystals put back, observations +8%, I/sigma 10.3 -> 10.6 and 10.7 -> 11.0, CC1/2 0.966 -> 0.969 and 0.970 -> 0.966. Tests: SublatticeRelabel_DoubledAxisGoesBackToTheTrueCell, _OtherSettingAndUnrelatedLattice, _SignalPicksTheAxisThatWasDoubled. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
866 lines
42 KiB
C++
866 lines
42 KiB
C++
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include <cstring>
|
|
#include <random>
|
|
#include "../image_analysis/bragg_integration/CalcISigma.h"
|
|
#include "../image_analysis/scale_merge/ErrorModel.h"
|
|
#include "../image_analysis/scale_merge/HKLKey.h"
|
|
#include "../image_analysis/scale_merge/Merge.h"
|
|
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
|
|
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
|
|
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
|
|
#include "../image_analysis/scale_merge/SublatticeRelabel.h"
|
|
#include "gemmi/reciproc.hpp"
|
|
|
|
TEST_CASE("HKLKey_NoSG_noMergeFriedel") {
|
|
HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1));
|
|
CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
|
|
}
|
|
|
|
TEST_CASE("HKLKey_NoSG_MergeFriedel") {
|
|
HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1));
|
|
CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
|
|
}
|
|
|
|
TEST_CASE("HKLKey_SG1_MergeFriedel") {
|
|
HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1));
|
|
CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
|
|
}
|
|
|
|
TEST_CASE("HKLKey_SG1_NoMergeFriedel") {
|
|
HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1));
|
|
CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
|
|
}
|
|
|
|
TEST_CASE("HKLKey_SG96_MergeFriedel") {
|
|
HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(96));
|
|
CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,1,3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,-2,3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,-1,3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(1,-2,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,2,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,1,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2, -1, -3));
|
|
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,-1,3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2, 1, 3));
|
|
}
|
|
|
|
TEST_CASE("HKLKey_SG96_NoMergeFriedel") {
|
|
HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(96));
|
|
CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3));
|
|
CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,1,3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,-2,3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,-1,3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(1,-2,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,2,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,1,-3));
|
|
CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2, -1, -3));
|
|
|
|
CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(-2,-1,3));
|
|
CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(2, 1, 3));
|
|
}
|
|
|
|
TEST_CASE("HKLKey_pack_friedel") {
|
|
HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1));
|
|
CHECK(hkl_key_gen(-1, -2, -3).pack() != hkl_key_gen(1,2,3).pack());
|
|
CHECK(hkl_key_gen(-1,-2,-3).pack() == hkl_key_gen(-1,-2,-3).pack());
|
|
CHECK(hkl_key_gen(-1,-2,-3).pack() != hkl_key_gen(1,-2,-3).pack());
|
|
}
|
|
|
|
TEST_CASE("HKLKey_pack_no_friedel") {
|
|
HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1));
|
|
CHECK(hkl_key_gen(-1, -2, -3).pack() == hkl_key_gen(1,2,3).pack());
|
|
CHECK(hkl_key_gen(-1,-2,-3).pack() == hkl_key_gen(-1,-2,-3).pack());
|
|
CHECK(hkl_key_gen(-1,-2,-3).pack() != hkl_key_gen(1,-2,-3).pack());
|
|
}
|
|
|
|
|
|
TEST_CASE("HKLKey_sys_absence_P212121") {
|
|
HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(19));
|
|
CHECK(hkl_key_gen.IsSystematicallyAbsent(5,0,0));
|
|
CHECK(!hkl_key_gen.IsSystematicallyAbsent(6,0,0));
|
|
CHECK(hkl_key_gen.IsSystematicallyAbsent(0,5,0));
|
|
CHECK(hkl_key_gen.IsSystematicallyAbsent(0,0,5));
|
|
CHECK(!hkl_key_gen.IsSystematicallyAbsent(0,4,0));
|
|
CHECK(!hkl_key_gen.IsSystematicallyAbsent(5,5,5));
|
|
}
|
|
TEST_CASE("AcceptReflection_ResolutionLimits") {
|
|
Reflection r{};
|
|
r.I = 100.0f;
|
|
r.sigma = 5.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.d = 20.0f;
|
|
|
|
// No limits: only the finiteness checks apply.
|
|
CHECK(AcceptReflection(r, std::nullopt, std::nullopt));
|
|
|
|
// Low-resolution limit rejects anything coarser than the limit, and is exclusive at it.
|
|
CHECK_FALSE(AcceptReflection(r, std::nullopt, std::optional<double>(15.0)));
|
|
CHECK(AcceptReflection(r, std::nullopt, std::optional<double>(20.0)));
|
|
CHECK(AcceptReflection(r, std::nullopt, std::optional<double>(50.0)));
|
|
|
|
// High-resolution limit still rejects anything finer, in the same direction as before.
|
|
CHECK_FALSE(AcceptReflection(r, std::optional<double>(25.0), std::nullopt));
|
|
CHECK(AcceptReflection(r, std::optional<double>(2.0), std::optional<double>(50.0)));
|
|
|
|
// The plain-double overload treats 0 as "no limit" at both ends.
|
|
CHECK(AcceptReflection(r, 0.0, 0.0));
|
|
CHECK_FALSE(AcceptReflection(r, 0.0, 15.0));
|
|
CHECK(AcceptReflection(r, 2.0, 50.0));
|
|
|
|
// A reflection with a saturated pixel is no measurement, whatever its resolution.
|
|
r.overloaded = true;
|
|
CHECK_FALSE(AcceptReflection(r, std::nullopt, std::nullopt));
|
|
CHECK_FALSE(AcceptReflection(r, 0.0, 0.0));
|
|
}
|
|
|
|
// --- Completeness denominator --------------------------------------------------------------------
|
|
|
|
namespace {
|
|
// A merged set built straight out of the reflections the cell and space group can give, so the
|
|
// test says exactly which of them were measured: every unique reflection between d_min and
|
|
// d_max_measured and none outside. Without Friedel merging an acentric contributes both hands.
|
|
std::vector<MergedReflection> MeasuredBetween(const gemmi::SpaceGroup &sg, const UnitCell &cell,
|
|
double d_min, double d_max_measured,
|
|
bool merge_friedel) {
|
|
const gemmi::UnitCell gemmi_cell = cell;
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
std::vector<MergedReflection> out;
|
|
for (const auto &hkl: gemmi::make_miller_vector(gemmi_cell, &sg, d_min, d_max_measured, true)) {
|
|
MergedReflection r;
|
|
r.h = hkl[0];
|
|
r.k = hkl[1];
|
|
r.l = hkl[2];
|
|
r.d = static_cast<float>(gemmi_cell.calculate_d(hkl));
|
|
r.I = 100.0f;
|
|
r.sigma = 10.0f;
|
|
r.I_half[0] = 100.0f;
|
|
r.I_half[1] = 100.0f;
|
|
out.push_back(r);
|
|
if (!merge_friedel && !gops.is_reflection_centric(hkl)) {
|
|
r.h = -hkl[0];
|
|
r.k = -hkl[1];
|
|
r.l = -hkl[2];
|
|
out.push_back(r);
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
|
|
MergeStatistics StatsWithLowLimit(const gemmi::SpaceGroup &sg, const std::optional<UnitCell> &cell,
|
|
const std::vector<MergedReflection> &merged,
|
|
std::optional<double> low_limit, bool merge_friedel) {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(sg);
|
|
ScalingSettings s = x.GetScalingSettings();
|
|
s.LowResolutionLimit_A(low_limit);
|
|
s.MergeFriedel(merge_friedel);
|
|
x.ImportScalingSettings(s);
|
|
|
|
MergeOnTheFly merge(x);
|
|
merge.ReferenceCell(cell);
|
|
return merge.MergeStats(merged, {});
|
|
}
|
|
|
|
double Completeness(const MergeStatisticsShell &s) {
|
|
return s.possible_unique_reflections > 0
|
|
? 100.0 * s.unique_reflections / s.possible_unique_reflections : 0.0;
|
|
}
|
|
|
|
const gemmi::SpaceGroup &TestSpaceGroup() { return gemmi::get_spacegroup_by_name("P 1 2 1"); }
|
|
constexpr UnitCell TEST_CELL{40, 45, 50, 90, 100, 90}; // synthetic; coarsest reflection ~49 A
|
|
}
|
|
|
|
// The low-resolution terms a beam stop ate must count as missing: the denominator is the declared
|
|
// range, so widening the declared range lowers completeness rather than leaving it alone.
|
|
TEST_CASE("MergeStats_CompletenessFallsWhenTheLowBoundCrossesAMaskedRegion") {
|
|
const auto &sg = TestSpaceGroup();
|
|
// Nothing coarser than 20 A was measured - it is all behind the stop.
|
|
const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);
|
|
REQUIRE(!merged.empty());
|
|
|
|
const auto at_20 = StatsWithLowLimit(sg, TEST_CELL, merged, 20.0, true);
|
|
const auto at_50 = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, true);
|
|
|
|
// Declared exactly where the data stop: everything possible was measured.
|
|
CHECK(Completeness(at_20.overall) > 99.0);
|
|
|
|
// Declared out to 50 A: the 20-50 A shell is in the denominator and in nothing else.
|
|
CHECK(at_50.overall.possible_unique_reflections > at_20.overall.possible_unique_reflections);
|
|
CHECK(at_50.overall.unique_reflections == at_20.overall.unique_reflections);
|
|
CHECK(Completeness(at_50.overall) < Completeness(at_20.overall));
|
|
// The innermost shell is where it bites.
|
|
CHECK(Completeness(at_50.shells.front()) < Completeness(at_20.shells.front()));
|
|
}
|
|
|
|
// No low-resolution limit means the whole sphere. The cell has no reflection coarser than 50 A, so
|
|
// freeing the 50 A limit must count the same set - the fix does not presuppose either default.
|
|
TEST_CASE("MergeStats_CompletenessWithNoDeclaredLowLimit") {
|
|
const auto &sg = TestSpaceGroup();
|
|
const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);
|
|
|
|
const auto at_50 = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, true);
|
|
const auto unlimited = StatsWithLowLimit(sg, TEST_CELL, merged, std::nullopt, true);
|
|
|
|
CHECK(unlimited.overall.possible_unique_reflections == at_50.overall.possible_unique_reflections);
|
|
CHECK(Completeness(unlimited.overall) == Catch::Approx(Completeness(at_50.overall)));
|
|
// The shell table stays finite even though the bound is not.
|
|
CHECK(std::isfinite(unlimited.shells.front().d_max));
|
|
}
|
|
|
|
// Counting the two Bijvoet mates of an acentric separately doubles the denominator too, so a fully
|
|
// measured anomalous set is 100% complete and not 200%.
|
|
TEST_CASE("MergeStats_CompletenessNeverExceeds100") {
|
|
const auto &sg = TestSpaceGroup();
|
|
for (const bool merge_friedel: {true, false}) {
|
|
const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 50.0, merge_friedel);
|
|
const auto stats = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, merge_friedel);
|
|
INFO("merge_friedel = " << merge_friedel);
|
|
CHECK(Completeness(stats.overall) <= 100.0);
|
|
CHECK(Completeness(stats.overall) > 99.0);
|
|
for (const auto &sh: stats.shells)
|
|
CHECK(Completeness(sh) <= 100.0);
|
|
}
|
|
}
|
|
|
|
// The shell grid and the denominator share their bounds, so every possible reflection lands in a
|
|
// shell: the shells sum to the overall, and the overall is the sphere the run declared.
|
|
TEST_CASE("MergeStats_PossibleSumsOverTheShellsToTheDeclaredSphere") {
|
|
const auto &sg = TestSpaceGroup();
|
|
const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);
|
|
const auto stats = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, true);
|
|
|
|
int sum = 0;
|
|
for (const auto &sh: stats.shells)
|
|
sum += sh.possible_unique_reflections;
|
|
CHECK(sum == stats.overall.possible_unique_reflections);
|
|
|
|
// Counted independently over the same declared range - nothing is lost between the two.
|
|
const gemmi::UnitCell gemmi_cell = TEST_CELL;
|
|
const int expected = gemmi::count_reflections(gemmi_cell, &sg, stats.overall.d_min * 0.999, 50.0, true);
|
|
CHECK(stats.overall.possible_unique_reflections == expected);
|
|
}
|
|
|
|
// Without a reference cell there is no set to count against; completeness stays unmeasured rather
|
|
// than becoming a number, with or without a declared low limit.
|
|
TEST_CASE("MergeStats_NoReferenceCellLeavesCompletenessUnmeasured") {
|
|
const auto &sg = TestSpaceGroup();
|
|
const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);
|
|
|
|
for (const std::optional<double> low_limit: {std::optional<double>(50.0), std::optional<double>()}) {
|
|
const auto stats = StatsWithLowLimit(sg, std::nullopt, merged, low_limit, true);
|
|
CHECK(stats.overall.possible_unique_reflections == 0);
|
|
CHECK(stats.overall.unique_reflections > 0);
|
|
CHECK(Completeness(stats.overall) == 0.0);
|
|
}
|
|
}
|
|
|
|
// The reference-range table (--report-resolution) is binned over the range it is GIVEN, whether or
|
|
// not the data reach it: its last shell ends at the declared d_min, and past the run's own limit
|
|
// nothing is read from anywhere - the reflections the run did not keep count as missing.
|
|
TEST_CASE("MergeStats_ReferenceRangeIsBinnedOverTheDeclaredRange") {
|
|
const auto &sg = TestSpaceGroup();
|
|
const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true); // what the run kept
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(sg);
|
|
MergeOnTheFly merge(x);
|
|
merge.ReferenceCell(TEST_CELL);
|
|
const auto own = merge.MergeStats(merged, {});
|
|
|
|
// A range coarser than the run's own: a complete subset, on a grid ending at the declared bound.
|
|
const auto within = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{3.0, 50.0});
|
|
CHECK(within.shells.back().d_min == Catch::Approx(3.0));
|
|
CHECK(within.shells.front().d_max == Catch::Approx(50.0));
|
|
int coarser_than_3 = 0;
|
|
for (const auto &m: merged)
|
|
if (m.d > 3.0f) ++coarser_than_3;
|
|
CHECK(within.overall.unique_reflections == coarser_than_3);
|
|
CHECK(within.overall.unique_reflections < own.overall.unique_reflections);
|
|
CHECK(Completeness(within.overall) > 99.0);
|
|
|
|
// A range finer than the run kept: the grid still ends at 1.5 A, the outer shell is empty and its
|
|
// possible reflections are in the denominator, and the measured range says where the data stop.
|
|
const auto beyond = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{1.5, 50.0});
|
|
CHECK(beyond.shells.back().d_min == Catch::Approx(1.5));
|
|
CHECK(beyond.overall.unique_reflections == own.overall.unique_reflections);
|
|
CHECK(beyond.overall.possible_unique_reflections > own.overall.possible_unique_reflections);
|
|
CHECK(Completeness(beyond.overall) < Completeness(own.overall));
|
|
CHECK(beyond.shells.back().unique_reflections == 0);
|
|
CHECK(beyond.shells.back().possible_unique_reflections > 0);
|
|
CHECK(beyond.overall.d_min == Catch::Approx(own.overall.d_min));
|
|
|
|
// A range the data never reach at all is still a grid, not an error.
|
|
const auto empty = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{1.0, 1.5});
|
|
CHECK(empty.overall.unique_reflections == 0);
|
|
CHECK(empty.overall.possible_unique_reflections > 0);
|
|
CHECK(empty.overall.d_min == 0.0f);
|
|
}
|
|
|
|
// ---------------------------------------------------------------- the automatic resolution cutoff
|
|
namespace {
|
|
// Half-set pairs spread uniformly in s = 1/d^2 over [s_from, s_to), either correlated with each
|
|
// other (signal) or drawn independently (noise, CC1/2 ~ 0), so a whole CC1/2 curve can be built
|
|
// band by band.
|
|
void AddBand(std::vector<MergedReflection> &v, std::mt19937 &rng,
|
|
double s_from, double s_to, int n, bool correlated, double sigma = 1.0) {
|
|
std::normal_distribution<double> g(0.0, sigma);
|
|
for (int j = 0; j < n; ++j) {
|
|
MergedReflection m;
|
|
m.h = static_cast<int32_t>(v.size()); // every reflection its own index
|
|
m.d = static_cast<float>(1.0 / std::sqrt(s_from + (j + 0.5) * (s_to - s_from) / n));
|
|
const double a = g(rng), b = g(rng);
|
|
m.I_half[0] = static_cast<float>(a);
|
|
m.I_half[1] = static_cast<float>(correlated ? a : b);
|
|
// The pair's weight in a CC1/2: the precision it would have had at the typical frame scale
|
|
// over the precision it has - 1/sigma^2 for a band scaled up by sigma from dead frames.
|
|
m.cc_weight = static_cast<float>(1.0 / (sigma * sigma));
|
|
v.push_back(m);
|
|
}
|
|
}
|
|
}
|
|
|
|
// A clean fall-off: CC1/2 crosses the target where the signal stops, and the cut is written one
|
|
// shell past it.
|
|
TEST_CASE("ResolutionCutoff_CleanFallOff") {
|
|
Logger logger("test");
|
|
std::mt19937 rng(12345);
|
|
std::vector<MergedReflection> merged;
|
|
AddBand(merged, rng, 0.01, 0.25, 480, true); // signal to 1/sqrt(0.25) = 2.00 A
|
|
AddBand(merged, rng, 0.25, 0.51, 520, false); // noise beyond it
|
|
|
|
const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
|
|
REQUIRE(rc.d_fit);
|
|
CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15));
|
|
REQUIRE(rc.d_cut);
|
|
CHECK(*rc.d_cut < *rc.d_fit); // the deliberate one-shell extension
|
|
CHECK(*rc.d_cut == Catch::Approx(1.92).margin(0.15));
|
|
}
|
|
|
|
// Past the fall-off, a few single reflections read hundreds of times the shell's noise in both halves -
|
|
// an artefact the half-sets reproduce. Each carries nearly all of its shell's covariance, so a Pearson
|
|
// CC1/2 reads those shells as signal and the cut runs on through them; one reflection per shell must not
|
|
// decide where the signal stops.
|
|
TEST_CASE("ResolutionCutoff_SingleReflectionsDoNotHoldTheCutOpen") {
|
|
Logger logger("test");
|
|
std::mt19937 rng(12345);
|
|
std::vector<MergedReflection> merged;
|
|
AddBand(merged, rng, 0.01, 0.25, 480, true);
|
|
AddBand(merged, rng, 0.25, 0.51, 520, false);
|
|
for (const double s : {0.26, 0.28, 0.30, 0.32}) {
|
|
MergedReflection m;
|
|
m.h = static_cast<int32_t>(merged.size());
|
|
m.d = static_cast<float>(1.0 / std::sqrt(s));
|
|
m.I_half[0] = m.I_half[1] = 300.0f;
|
|
merged.push_back(m);
|
|
}
|
|
|
|
const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
|
|
REQUIRE(rc.d_fit);
|
|
CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15));
|
|
}
|
|
|
|
// A sweep with a long stretch where the crystal barely diffracted: the reflections measured only there
|
|
// are noise scaled up by 1/G, sigmas with them, and there are many of them at every resolution. Counted
|
|
// by the information they carry they must not hide where the well-measured reflections stop.
|
|
TEST_CASE("ResolutionCutoff_ScaledUpNoiseDoesNotHideTheFallOff") {
|
|
Logger logger("test");
|
|
std::mt19937 rng(12345);
|
|
std::vector<MergedReflection> merged;
|
|
AddBand(merged, rng, 0.01, 0.25, 480, true);
|
|
AddBand(merged, rng, 0.25, 0.51, 520, false);
|
|
AddBand(merged, rng, 0.01, 0.51, 300, false, 50.0);
|
|
|
|
const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
|
|
REQUIRE(rc.d_fit);
|
|
CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15));
|
|
|
|
// Counted as equals, the same reflections hide it: the curve reads noise from the first bin on.
|
|
for (auto &m : merged) m.cc_weight = 1.0f;
|
|
const auto unweighted = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
|
|
CHECK((!unweighted.d_fit || *unweighted.d_fit > 3.0));
|
|
}
|
|
|
|
// The frame factor of the CC1/2 weight is 1 on a sweep whose frames all sit at one scale, however that
|
|
// scale is spread over the frames, and 1 for any frame brighter than typical.
|
|
TEST_CASE("CCHalfFrameFactors_UniformScaleIsOne") {
|
|
const auto f = CCHalfFrameFactors({1.7, 1.7, 1.7, 1.7}, {100, 3, 250, 0});
|
|
for (double x : f) CHECK(x == 1.0);
|
|
const auto g = CCHalfFrameFactors({1.0, 2.0}, {100, 100});
|
|
CHECK(g[1] == 1.0);
|
|
CHECK(g[0] == Catch::Approx(9.0 / 5.0 * 9.0 / 5.0)); // G_ref = (1 + 8) / (1 + 4)
|
|
}
|
|
|
|
// Most of the sweep at 2% of the good frames' scale: the typical scale is still the good frames', not
|
|
// the run median, and an observation from a dead frame carries 2500x the variance of a good one.
|
|
TEST_CASE("CCHalfFrameFactors_DeadStretchDoesNotSetTheTypicalScale") {
|
|
std::vector<double> scale;
|
|
std::vector<int64_t> n;
|
|
for (int f = 0; f < 100; ++f) {
|
|
scale.push_back(f < 60 ? 0.02 : 1.0);
|
|
n.push_back(50);
|
|
}
|
|
scale.push_back(NAN); // a frame without a scale counts as G = 1
|
|
n.push_back(50);
|
|
const auto factor = CCHalfFrameFactors(scale, n);
|
|
CHECK(factor[0] == Catch::Approx(2500.0).epsilon(0.01));
|
|
CHECK(factor[99] == Catch::Approx(1.0).epsilon(1e-3));
|
|
CHECK(factor[100] == Catch::Approx(1.0).epsilon(1e-3));
|
|
}
|
|
|
|
// Equal weights give the ordinary R_meas terms, whatever their common value.
|
|
TEST_CASE("WeightedRmeas_EqualWeightsAreTheOrdinaryRmeas") {
|
|
// Three observations 9, 10, 11 of a reflection with <I> = 10: sum|dev| = 2, sum I = 30.
|
|
double num, den;
|
|
REQUIRE(WeightedRmeasTerms(0.5 * 2.0, 0.5 * 30.0, 0.5 * 3, 0.25 * 3, num, den));
|
|
CHECK(num == Catch::Approx(std::sqrt(3.0 / 2.0) * 2.0));
|
|
CHECK(den == Catch::Approx(30.0));
|
|
}
|
|
|
|
// An observation that carries almost no information counts for almost nothing: the terms tend to those
|
|
// of the reflection without it, and a reflection left with one effective observation has none.
|
|
TEST_CASE("WeightedRmeas_NegligibleWeightDropsOut") {
|
|
// 9 and 11 at weight 1, <I> = 10; a third at 40 with weight 1e-6.
|
|
const double v = 1e-6;
|
|
double num, den;
|
|
REQUIRE(WeightedRmeasTerms(2.0 + v * 30.0, 20.0 + v * 40.0, 2.0 + v, 2.0 + v * v, num, den));
|
|
CHECK(num == Catch::Approx(std::sqrt(2.0) * 2.0).epsilon(1e-4));
|
|
CHECK(den == Catch::Approx(20.0).epsilon(1e-4));
|
|
CHECK_FALSE(WeightedRmeasTerms(0.0, 10.0, 1.0, 1.0, num, den));
|
|
}
|
|
|
|
// A fall-off region a logistic cannot follow: CC1/2 drops through the target and comes straight back
|
|
// up. The fitted crossing is then an extrapolation far past the bins it was made over, and reading
|
|
// the cut off it writes the data deep into the noise; the crossing the bins themselves show is where
|
|
// the signal stopped, and that is what must be used.
|
|
TEST_CASE("ResolutionCutoff_RaggedFallOffIsReadOffTheBins") {
|
|
Logger logger("test");
|
|
std::mt19937 rng(12345);
|
|
std::vector<MergedReflection> merged;
|
|
AddBand(merged, rng, 0.01, 0.13, 240, true); // signal to 1/sqrt(0.13) = 2.77 A
|
|
AddBand(merged, rng, 0.13, 0.17, 80, false); // a hole below the target
|
|
AddBand(merged, rng, 0.17, 0.25, 160, true); // correlated again - not a fall-off
|
|
AddBand(merged, rng, 0.25, 0.51, 520, false);
|
|
|
|
const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
|
|
REQUIRE(rc.d_fit);
|
|
CHECK(*rc.d_fit == Catch::Approx(2.77).margin(0.20));
|
|
REQUIRE(rc.d_cut);
|
|
CHECK(*rc.d_cut > 2.30); // coarser than the band that correlates again
|
|
}
|
|
|
|
namespace {
|
|
// Samples of var = a*s2 + b^2*I2 over four decades of counting I/sigma, with counting variances
|
|
// spread over a factor of 100 at every intensity, and a fraction of gross outliers.
|
|
std::vector<ErrorModelSample> SyntheticErrorModelSamples(double a, double b, double max_snr,
|
|
double outlier_fraction) {
|
|
std::mt19937 rng(7);
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
std::normal_distribution<double> z(0.0, 1.0);
|
|
std::vector<ErrorModelSample> out;
|
|
for (int i = 0; i < 200000; ++i) {
|
|
const double s2 = std::pow(10.0, 2.0 * u(rng));
|
|
const double snr = max_snr * std::pow(10.0, -4.0 * u(rng));
|
|
const double I2 = snr * snr * s2;
|
|
const double e = z(rng);
|
|
double dev2 = (a * s2 + b * b * I2) * e * e;
|
|
if (u(rng) < outlier_fraction) dev2 *= 400.0;
|
|
out.push_back({s2, I2, dev2, 2.0f});
|
|
}
|
|
return out;
|
|
}
|
|
}
|
|
|
|
// The fit recovers a and b from observations whose counting variances differ widely inside every
|
|
// intensity bin (where a ratio of bin medians does not), with 0.3% gross outliers in the pool.
|
|
TEST_CASE("ErrorModel_RecoversAAndBThroughOutliers") {
|
|
const auto pool = SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.003);
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto fit = FitErrorModel(pool, scratch, 4);
|
|
REQUIRE(fit.active);
|
|
CHECK(fit.b_measured);
|
|
CHECK(fit.b_resolved);
|
|
CHECK(fit.a == Catch::Approx(1.3).epsilon(0.03));
|
|
CHECK(1.0 / std::sqrt(fit.b2) == Catch::Approx(1.0 / 0.03).epsilon(0.05));
|
|
}
|
|
|
|
// A few observations whose counting variance is hugely overstated - their scatter is ordinary - do not
|
|
// set the fit: each sample counts by its deviation relative to its own variance.
|
|
TEST_CASE("ErrorModel_OverstatedCountingVarianceDoesNotSetTheFit") {
|
|
auto pool = SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.0);
|
|
for (size_t i = 0; i < pool.size(); i += 200)
|
|
pool[i].s2 *= 1e4;
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto fit = FitErrorModel(pool, scratch, 4);
|
|
REQUIRE(fit.active);
|
|
CHECK(fit.a == Catch::Approx(1.3).epsilon(0.03));
|
|
CHECK(1.0 / std::sqrt(fit.b2) == Catch::Approx(1.0 / 0.03).epsilon(0.05));
|
|
}
|
|
|
|
// The same pool on one thread and on many gives the same numbers.
|
|
TEST_CASE("ErrorModel_SameOnAnyNumberOfThreads") {
|
|
const auto pool = SyntheticErrorModelSamples(1.1, 0.05, 100.0, 0.001);
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto one = FitErrorModel(pool, scratch, 1);
|
|
const auto many = FitErrorModel(pool, scratch, 16);
|
|
CHECK(one.a == many.a);
|
|
CHECK(one.b2 == many.b2);
|
|
}
|
|
|
|
// Data that never reach where b could be seen: a is fitted alone and b held at 0.
|
|
TEST_CASE("ErrorModel_BNotMeasuredOnWeakData") {
|
|
const auto pool = SyntheticErrorModelSamples(0.9, 0.05, 1.5, 0.0);
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto fit = FitErrorModel(pool, scratch, 4);
|
|
REQUIRE(fit.active);
|
|
CHECK(!fit.b_measured);
|
|
CHECK(fit.b2 == 0.0);
|
|
CHECK(fit.a == Catch::Approx(0.9).epsilon(0.03));
|
|
}
|
|
|
|
#ifdef JFJOCH_USE_CUDA
|
|
#include "../common/CUDAWrapper.h"
|
|
#include "../image_analysis/scale_merge/ErrorModelGPU.h"
|
|
|
|
// The device fit against the host one. Same bins and same medians; only the order each bin's two sums
|
|
// are added in differs, so a and b^2 agree to rounding and every flag agrees exactly. The same pool
|
|
// twice on the device gives the same bits.
|
|
TEST_CASE("ErrorModel_DeviceFitIsTheHostFit") {
|
|
if (get_gpu_count() == 0) {
|
|
WARN("No CUDA GPU present. Skipping ErrorModel_DeviceFitIsTheHostFit");
|
|
return;
|
|
}
|
|
for (const auto &pool : {SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.003),
|
|
SyntheticErrorModelSamples(0.9, 0.05, 1.5, 0.0)}) {
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto host = FitErrorModel(pool, scratch, 4);
|
|
const auto dev = FitErrorModelGPU(pool);
|
|
const auto again = FitErrorModelGPU(pool);
|
|
CHECK(dev.active == host.active);
|
|
CHECK(dev.b_measured == host.b_measured);
|
|
CHECK(dev.b_resolved == host.b_resolved);
|
|
CHECK(dev.a == Catch::Approx(host.a).epsilon(1e-12));
|
|
CHECK(dev.b2 == Catch::Approx(host.b2).epsilon(1e-12).margin(1e-300));
|
|
CHECK(again.a == dev.a);
|
|
CHECK(again.b2 == dev.b2);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
// The dataset-wide Wilson B: a pure exp(-B / 2d^2) fall-off, read back over the fit window. A crystal
|
|
// diffracting to just under 4 A has almost nothing between its d_min and the 4 A convention, so the
|
|
// window is widened to the 1/d^2 width a 3 A dataset gets; and a high B is a measurement, not garbage.
|
|
TEST_CASE("Global Wilson B over a widened window, with no upper cap", "[wilson]") {
|
|
const auto merge = [](double b, double d_min) {
|
|
std::vector<MergedReflection> out;
|
|
for (int i = 0; i < 20000; ++i) {
|
|
// Uniform in 1/d^3 from 20 A to d_min, as reflections are.
|
|
const double s3 = (i + 0.5) / 20000.0 * (1.0 / std::pow(d_min, 3) - 1.0 / 8000.0) + 1.0 / 8000.0;
|
|
MergedReflection r;
|
|
r.h = i; r.k = 0; r.l = 0;
|
|
r.d = static_cast<float>(1.0 / std::cbrt(s3));
|
|
r.I = static_cast<float>(1e6 * std::exp(-b / (2.0 * r.d * r.d)));
|
|
r.sigma = 0.01f * r.I;
|
|
out.push_back(r);
|
|
}
|
|
return out;
|
|
};
|
|
const GlobalWilsonB high_res = CalcGlobalWilsonB(merge(40.0, 2.0));
|
|
CHECK(high_res.b == Catch::Approx(40.0).epsilon(0.02));
|
|
CHECK(high_res.d_low_A == Catch::Approx(4.0));
|
|
|
|
const GlobalWilsonB low_res = CalcGlobalWilsonB(merge(300.0, 3.9));
|
|
CHECK(low_res.b == Catch::Approx(300.0).epsilon(0.02));
|
|
CHECK(low_res.n_shells >= 3);
|
|
CHECK(1.0 / (3.9 * 3.9) - 1.0 / (low_res.d_low_A * low_res.d_low_A)
|
|
== Catch::Approx(1.0 / 9.0 - 1.0 / 16.0).epsilon(0.01));
|
|
}
|
|
|
|
// The parallel AddImages is the serial AddImage loop, bit for bit: the per-observation work is spread
|
|
// over the workers, but every sum and every half-set choice is still made in image order. More images
|
|
// than one of its blocks, so the block boundary is crossed.
|
|
TEST_CASE("MergeOnTheFly_AddImagesIsTheSerialMerge") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(TestSpaceGroup());
|
|
|
|
std::mt19937 rng(7);
|
|
std::uniform_int_distribution<int> index(-6, 6);
|
|
std::uniform_real_distribution<float> unit(0.0f, 1.0f);
|
|
std::vector<IntegrationOutcome> outcomes(5000);
|
|
for (auto &o: outcomes) {
|
|
for (int j = 0; j < 20; ++j) {
|
|
Reflection r{};
|
|
r.h = index(rng);
|
|
r.k = index(rng);
|
|
r.l = index(rng);
|
|
if (r.h == 0 && r.k == 0 && r.l == 0)
|
|
continue;
|
|
r.d = 20.0f / std::sqrt(static_cast<float>(r.h * r.h + r.k * r.k + r.l * r.l));
|
|
r.I = 100.0f * unit(rng) - 5.0f;
|
|
r.sigma = 1.0f + 5.0f * unit(rng);
|
|
r.var_bkg = r.sigma * r.sigma * unit(rng);
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 0.5f + 0.5f * unit(rng);
|
|
r.image_scale_corr = 0.5f + unit(rng);
|
|
r.on_ice_ring = unit(rng) < 0.05f;
|
|
o.reflections.push_back(r);
|
|
}
|
|
}
|
|
|
|
MergeOnTheFly serial(x);
|
|
serial.RefineErrorModel(outcomes);
|
|
for (size_t i = 0; i < outcomes.size(); ++i)
|
|
serial.AddImage(outcomes[i], static_cast<int64_t>(i));
|
|
const auto serial_merged = serial.ExportReflections();
|
|
|
|
MergeOnTheFly parallel(x);
|
|
parallel.Threads(8);
|
|
parallel.RefineErrorModel(outcomes);
|
|
parallel.AddImages(outcomes);
|
|
const auto parallel_merged = parallel.ExportReflections();
|
|
|
|
REQUIRE(serial.ErrorModelActive());
|
|
CHECK(serial.ErrorModelA() == parallel.ErrorModelA());
|
|
CHECK(serial.ErrorModelB() == parallel.ErrorModelB());
|
|
REQUIRE(serial_merged.size() == parallel_merged.size());
|
|
for (size_t i = 0; i < serial_merged.size(); ++i) {
|
|
CHECK(serial_merged[i].h == parallel_merged[i].h);
|
|
CHECK(serial_merged[i].k == parallel_merged[i].k);
|
|
CHECK(serial_merged[i].l == parallel_merged[i].l);
|
|
CHECK(serial_merged[i].I == parallel_merged[i].I);
|
|
CHECK(serial_merged[i].sigma == parallel_merged[i].sigma);
|
|
for (int half = 0; half < 2; ++half) // NaN where a half is empty, on both sides
|
|
CHECK(std::memcmp(&serial_merged[i].I_half[half], ¶llel_merged[i].I_half[half], sizeof(float)) == 0);
|
|
}
|
|
|
|
const auto serial_stats = serial.MergeStats(serial_merged, outcomes);
|
|
const auto parallel_stats = parallel.MergeStats(parallel_merged, outcomes);
|
|
CHECK(serial_stats.overall.total_observations == parallel_stats.overall.total_observations);
|
|
CHECK(serial_stats.overall.cc_half == parallel_stats.overall.cc_half);
|
|
CHECK(serial_stats.overall.r_meas == parallel_stats.overall.r_meas);
|
|
}
|
|
|
|
// Stills partiality post-refinement on noiseless synthetic crystals: each crystal's true orientation is
|
|
// tilted from the stored one by up to 0.05 deg, its intensities are the true partials of a fixed set of
|
|
// full intensities, and its stored partialities those of the stored orientation. The refinement has to
|
|
// move the partialities towards the true ones. (The Ceres solve this replaced gave the same numbers: error
|
|
// 0.0559 -> 0.0302, mean tilt 0.0206 deg.)
|
|
TEST_CASE("StillsPartialityRefine_RecoversTheTilt") {
|
|
constexpr float CELL = 50.0f, WAVELENGTH = 1.0f, D_MIN = 2.0f;
|
|
constexpr double GAMMA0 = 0.001, ACCEPT = 1.7 * GAMMA0; // RLP width; |dist_ewald| accepted at prediction
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
DiffractionGeometry geom;
|
|
geom.Wavelength_A(WAVELENGTH);
|
|
const Coord S0 = geom.GetScatteringVector();
|
|
|
|
std::mt19937 rng(11);
|
|
std::normal_distribution<double> gauss(0.0, 1.0);
|
|
std::uniform_real_distribution<double> unit(0.0, 1.0);
|
|
auto full_I = [](int h, int k, int l) { // a fixed, Friedel-symmetric pseudo-random full intensity
|
|
const uint64_t key = HKLKeyGenerator(true, gemmi::get_spacegroup_by_name("P 1"))(h, k, l).pack();
|
|
return 10.0 + 990.0 * static_cast<double>((key * 0x9e3779b97f4a7c15ULL) >> 40) / static_cast<double>(1ULL << 24);
|
|
};
|
|
auto rotate = [](const double R[9], const Coord &v) {
|
|
return Coord(static_cast<float>(R[0] * v.x + R[1] * v.y + R[2] * v.z),
|
|
static_cast<float>(R[3] * v.x + R[4] * v.y + R[5] * v.z),
|
|
static_cast<float>(R[6] * v.x + R[7] * v.y + R[8] * v.z));
|
|
};
|
|
auto rotation = [](double ax, double ay, double az, double angle, double R[9]) { // Rodrigues
|
|
const double n = std::sqrt(ax * ax + ay * ay + az * az);
|
|
ax /= n; ay /= n; az /= n;
|
|
const double c = std::cos(angle), s = std::sin(angle), t = 1.0 - c;
|
|
const double M[9] = {c + t * ax * ax, t * ax * ay - s * az, t * ax * az + s * ay,
|
|
t * ax * ay + s * az, c + t * ay * ay, t * ay * az - s * ax,
|
|
t * ax * az - s * ay, t * ay * az + s * ax, c + t * az * az};
|
|
std::copy(M, M + 9, R);
|
|
};
|
|
auto partiality = [&](const Coord &q) {
|
|
const double de = (q + S0).Length() - 1.0 / WAVELENGTH;
|
|
return std::exp(-0.5 * de * de / (GAMMA0 * GAMMA0));
|
|
};
|
|
|
|
std::vector<IntegrationOutcome> outcomes(150);
|
|
std::vector<std::vector<double>> true_p(outcomes.size());
|
|
const int hmax = static_cast<int>(CELL / D_MIN);
|
|
for (size_t c = 0; c < outcomes.size(); ++c) {
|
|
double R[9], T[9];
|
|
rotation(gauss(rng), gauss(rng), gauss(rng), 2.0 * M_PI * unit(rng), R);
|
|
const double tilt = (0.02 + 0.03 * unit(rng)) * M_PI / 180.0, phi = 2.0 * M_PI * unit(rng);
|
|
rotation(std::cos(phi), std::sin(phi), 0.0, tilt, T); // about an axis normal to the beam
|
|
auto &o = outcomes[c];
|
|
o.geom = geom;
|
|
o.latt = CrystalLattice(rotate(R, Coord(CELL, 0, 0)), rotate(R, Coord(0, CELL, 0)),
|
|
rotate(R, Coord(0, 0, CELL)));
|
|
const double G = 0.5 + unit(rng);
|
|
for (int h = -hmax; h <= hmax; ++h)
|
|
for (int k = -hmax; k <= hmax; ++k)
|
|
for (int l = -hmax; l <= hmax; ++l) {
|
|
const Coord q = o.latt.Astar() * h + o.latt.Bstar() * k + o.latt.Cstar() * l;
|
|
const float dstar = q.Length();
|
|
if (dstar == 0.0f || dstar > 1.0f / D_MIN)
|
|
continue;
|
|
if (std::fabs((q + S0).Length() - 1.0 / WAVELENGTH) > ACCEPT)
|
|
continue;
|
|
const double p_true = partiality(rotate(T, q));
|
|
Reflection r{};
|
|
r.h = h; r.k = k; r.l = l;
|
|
r.d = 1.0f / dstar;
|
|
r.I = static_cast<float>(G * p_true * full_I(h, k, l));
|
|
r.sigma = static_cast<float>(1.0 + 0.01 * r.I);
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = static_cast<float>(partiality(q));
|
|
r.image_scale_corr = 1.0f / r.partiality;
|
|
o.reflections.push_back(r);
|
|
true_p[c].push_back(p_true);
|
|
}
|
|
REQUIRE(o.reflections.size() >= 40);
|
|
}
|
|
|
|
auto mean_abs_error = [&]() {
|
|
double sum = 0.0;
|
|
size_t n = 0;
|
|
for (size_t c = 0; c < outcomes.size(); ++c)
|
|
for (size_t i = 0; i < outcomes[c].reflections.size(); ++i, ++n)
|
|
sum += std::fabs(outcomes[c].reflections[i].partiality - true_p[c][i]);
|
|
return sum / static_cast<double>(n);
|
|
};
|
|
const double before = mean_abs_error();
|
|
ScaleOnTheFly(x, MergeAll(x, outcomes, 4, std::nullopt)).Scale(outcomes, 4);
|
|
const double mean_tilt = StillsPartialityRefine(x).Run(outcomes, 4, std::nullopt);
|
|
const double after = mean_abs_error();
|
|
CHECK(mean_tilt > 0.015);
|
|
CHECK(mean_tilt < 0.06);
|
|
CHECK(after < 0.6 * before); // measured 0.54: the prior holds the larger tilts back
|
|
}
|
|
|
|
namespace {
|
|
// The basis vectors a, b, c of a cell, turned by a fixed arbitrary rotation, as a frame's lattice.
|
|
std::array<Coord, 3> TurnedAxes(const UnitCell &cell) {
|
|
const CrystalLattice l(cell);
|
|
const double c = std::cos(0.7), s = std::sin(0.7); // about (1, 2, 3)/sqrt(14)
|
|
const double n[3] = {1.0 / std::sqrt(14.0), 2.0 / std::sqrt(14.0), 3.0 / std::sqrt(14.0)};
|
|
auto turn = [&](const Coord &v) {
|
|
const double dot = n[0] * v.x + n[1] * v.y + n[2] * v.z;
|
|
const double x[3] = {n[1] * v.z - n[2] * v.y, n[2] * v.x - n[0] * v.z, n[0] * v.y - n[1] * v.x};
|
|
return Coord(static_cast<float>(v.x * c + x[0] * s + n[0] * dot * (1 - c)),
|
|
static_cast<float>(v.y * c + x[1] * s + n[1] * dot * (1 - c)),
|
|
static_cast<float>(v.z * c + x[2] * s + n[2] * dot * (1 - c)));
|
|
};
|
|
return {turn(l.Vec0()), turn(l.Vec1()), turn(l.Vec2())};
|
|
}
|
|
|
|
Coord ScatteringVector(const CrystalLattice &l, const Reflection &r) {
|
|
return l.Astar() * static_cast<float>(r.h) + l.Bstar() * static_cast<float>(r.k)
|
|
+ l.Cstar() * static_cast<float>(r.l);
|
|
}
|
|
|
|
IntegrationOutcome FrameWithAllIndices(const CrystalLattice &latt, int hmax) {
|
|
IntegrationOutcome o;
|
|
o.latt = latt;
|
|
for (int h = -hmax; h <= hmax; ++h)
|
|
for (int k = -hmax; k <= hmax; ++k)
|
|
for (int l = -hmax; l <= hmax; ++l) {
|
|
Reflection r{};
|
|
r.h = h; r.k = k; r.l = l;
|
|
r.predicted_x = static_cast<float>(o.reflections.size()); // identifies the observation
|
|
r.I = 100.0f;
|
|
r.sigma = 10.0f;
|
|
const float q = ScatteringVector(latt, r).Length();
|
|
r.d = q > 0.0f ? 1.0f / q : 100.0f;
|
|
o.reflections.push_back(r);
|
|
}
|
|
return o;
|
|
}
|
|
}
|
|
|
|
// A frame indexed with its c axis doubled: its reciprocal lattice holds every node of the true one plus
|
|
// the planes in between. Relabelling keeps exactly the true nodes, each at the same scattering vector,
|
|
// and gives the frame the true cell.
|
|
TEST_CASE("SublatticeRelabel_DoubledAxisGoesBackToTheTrueCell") {
|
|
const UnitCell reference{60, 60, 25, 90, 90, 90}; // synthetic tetragonal
|
|
const auto ax = TurnedAxes(reference);
|
|
const CrystalLattice super(ax[0], ax[1], ax[2] * 2.0f);
|
|
// A frame on the true cell beside it sets how fine the relabelled reflections may go.
|
|
std::vector<IntegrationOutcome> outcomes{FrameWithAllIndices(super, 4),
|
|
FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2]), 4)};
|
|
const auto before = outcomes[0];
|
|
|
|
IndexingSettings settings;
|
|
const auto result = RelabelOnSublattice(outcomes, reference, settings);
|
|
CHECK(result.frames == 1);
|
|
CHECK(result.reflections_kept == 9 * 9 * 5); // l even: -4, -2, 0, 2, 4
|
|
CHECK(result.reflections_dropped == 9 * 9 * 4);
|
|
CHECK(outcomes[0].latt.GetUnitCell().is_close(reference, 0.001, 0.01));
|
|
|
|
for (const auto &r: outcomes[0].reflections) {
|
|
const auto it = std::find_if(before.reflections.begin(), before.reflections.end(),
|
|
[&](const Reflection &b) { return b.predicted_x == r.predicted_x; });
|
|
REQUIRE(it != before.reflections.end());
|
|
CHECK(it->l % 2 == 0);
|
|
const Coord q_before = ScatteringVector(before.latt, *it);
|
|
const Coord q_after = ScatteringVector(outcomes[0].latt, r);
|
|
CHECK((q_after - q_before).Length() < 1e-5f);
|
|
}
|
|
}
|
|
|
|
// The true lattice in another setting (axes permuted) is the same lattice: every reflection is kept and
|
|
// only re-indexed. A frame whose lattice does not contain the reference is left alone.
|
|
TEST_CASE("SublatticeRelabel_OtherSettingAndUnrelatedLattice") {
|
|
const UnitCell reference{40, 50, 60, 90, 90, 90}; // synthetic orthorhombic
|
|
const auto ax = TurnedAxes(reference);
|
|
std::vector<IntegrationOutcome> outcomes{FrameWithAllIndices(CrystalLattice(ax[1], ax[2], ax[0]), 3),
|
|
FrameWithAllIndices(CrystalLattice(UnitCell{37, 53, 71, 90, 90, 90}), 3),
|
|
FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2]), 4)};
|
|
const auto permuted = outcomes[0];
|
|
const auto unrelated = outcomes[1];
|
|
|
|
IndexingSettings settings;
|
|
const auto result = RelabelOnSublattice(outcomes, reference, settings);
|
|
CHECK(result.frames == 1);
|
|
CHECK(result.reflections_dropped == 0);
|
|
CHECK(outcomes[0].latt.GetUnitCell().is_close(reference, 0.001, 0.01));
|
|
for (size_t i = 0; i < outcomes[0].reflections.size(); ++i) {
|
|
const Coord q_before = ScatteringVector(permuted.latt, permuted.reflections[i]);
|
|
const Coord q_after = ScatteringVector(outcomes[0].latt, outcomes[0].reflections[i]);
|
|
CHECK((q_after - q_before).Length() < 1e-5f);
|
|
}
|
|
CHECK(outcomes[1].reflections.size() == unrelated.reflections.size());
|
|
CHECK(outcomes[1].reflections[5].h == unrelated.reflections[5].h);
|
|
}
|
|
|
|
// A doubled axis that makes the supercell metric cubic fits the reference along any of the three axes;
|
|
// only the true one has signal at its nodes, and that is the one taken.
|
|
TEST_CASE("SublatticeRelabel_SignalPicksTheAxisThatWasDoubled") {
|
|
const UnitCell reference{60, 60, 30, 90, 90, 90}; // synthetic tetragonal; doubled c gives a cube
|
|
const auto ax = TurnedAxes(reference);
|
|
std::vector<IntegrationOutcome> outcomes{FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2] * 2.0f), 4),
|
|
FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2]), 4)};
|
|
for (auto &r: outcomes[0].reflections)
|
|
r.I = r.l % 2 == 0 ? 100.0f : 0.0f; // the nodes between the true ones are empty
|
|
const auto before = outcomes[0];
|
|
|
|
IndexingSettings settings;
|
|
const auto result = RelabelOnSublattice(outcomes, reference, settings);
|
|
REQUIRE(result.frames == 1);
|
|
for (const auto &r: outcomes[0].reflections) {
|
|
const auto &b = before.reflections[static_cast<size_t>(r.predicted_x)];
|
|
CHECK(b.l % 2 == 0);
|
|
CHECK((ScatteringVector(outcomes[0].latt, r) - ScatteringVector(before.latt, b)).Length() < 1e-5f);
|
|
}
|
|
}
|