leonarski_f
ced02f1f78
A frame whose background region reads below zero - a detector module whose pedestal is off on that
frame - integrates reflections with a negative background variance, a sigma floored at one photon and a
net intensity of thousands (I/sigma ~8000). The unscaled merge the per-image scales are fitted against
weights such an observation ~100x above any real one, so it takes its reflection over; every frame that
sees that reflection and finds nothing there is then pulled towards a scale of 0 (the proportional model
error caps such an observation's normalised residual at ~3.3 sigma, so the Cauchy fit never rejects it
and its pull grows as G falls), lands under the 1/50 floor and is dropped, and the partiality refit
against the survivors drops more.
On a serial run whose detector shows this on one module every 324th frame, always on the same residue
of the frame number modulo 4, two same-arm quarters of the run (frames mod 4 = 0 and = 2) merged to
CC1/2 0.78 and 0.28. The collapsed frames were good frames (CC with the deposited model's Fcalc 0.29
against 0.34 for the kept ones, the same hand) whose CC with the scaling reference was 0.03: the
reference, not the frames, was wrong. Measured (frames mod 4 = 0 / = 2, run's geometry):
before: 446 / 1355 frames dropped, R_split 32.6 / 64.0%, <I^2>/<I>^2 2.08 / 2.74, <|L|> 0.475 / 0.788,
CC of the merge with the model Fcalc (d >= 2.5 A) 0.842 / 0.381
after: 446 / 412 dropped, R_split 32.6 / 35.9%, 2.08 / 2.17, 0.475 / 0.478, 0.842 / 0.816
The whole run: 9081 -> 3081 frames dropped, CC1/2 0.59 -> 0.955, R_split 18.6 -> 14.9%, Fcalc CC
0.879 -> 0.884 (d >= 2.5 A); two runs without the defect are unchanged.
AcceptReflection now refuses such an observation (a NaN variance, from older files, is kept), so it
leaves every stills scale fit, merge and ambiguity correlation, and the stills merge logs how many
observations on how many frames were left out.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
(cherry picked from commit bd5c59915fd7698afcc0ef81722e562c78c5ee3a)
1641 lines
78 KiB
C++
1641 lines
78 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));
|
|
|
|
// Nor is one whose background variance is below zero; an unknown one (NaN, older files) is kept.
|
|
r.overloaded = false;
|
|
r.var_bkg = -1.0f;
|
|
CHECK_FALSE(AcceptReflection(r, std::nullopt, std::nullopt));
|
|
CHECK_FALSE(AcceptReflection(r, 0.0, 0.0));
|
|
r.var_bkg = NAN;
|
|
CHECK(AcceptReflection(r, std::nullopt, std::nullopt));
|
|
}
|
|
|
|
// --- 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;
|
|
}
|
|
|
|
// The statistics over a declared low-resolution bound (--report-resolution), or over the whole sphere
|
|
// without one: the scaling's own low-resolution limit bounds the scale fit, not the merge.
|
|
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.MergeFriedel(merge_friedel);
|
|
x.ImportScalingSettings(s);
|
|
|
|
MergeOnTheFly merge(x);
|
|
merge.ReferenceCell(cell);
|
|
const auto own = merge.MergeStats(merged, {});
|
|
if (!low_limit)
|
|
return own;
|
|
// Down to the data's own d_min, padded as MergeStats pads it.
|
|
return merge.MergeStats(merged, {}, {}, std::nullopt,
|
|
ReportResolutionRange{.d_min = own.overall.d_min * 0.999, .d_max = *low_limit});
|
|
}
|
|
|
|
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.
|
|
// The plain fit (the quoted FITTED_RESOLUTION) is that, with the weights left in place.
|
|
const auto plain = ComputeCCHalfLogisticCutoff(merged, 0.30, logger, /*weighted=*/false);
|
|
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));
|
|
CHECK(plain.d_fit == unweighted.d_fit);
|
|
CHECK(plain.d_cut == unweighted.d_cut);
|
|
}
|
|
|
|
// 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, i / 8}); // eight mates per reflection
|
|
}
|
|
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);
|
|
}
|
|
|
|
// The stills model: samples drawn with a known a, b and partiality term c - c^2 ((1 - p) / p)^2 <I>^2, p spread
|
|
// over [0.1, 1] - give those three back.
|
|
TEST_CASE("StillsErrorModel_RecoversThePartialityTerm") {
|
|
std::mt19937 rng(7);
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
std::normal_distribution<double> z(0.0, 1.0);
|
|
const double a = 1.2, b = 0.08, c = 0.15;
|
|
std::vector<StillsErrorModelSample> pool;
|
|
for (int i = 0; i < 200000; ++i) {
|
|
const double s2 = std::pow(10.0, 2.0 * u(rng));
|
|
const double snr = 300.0 * std::pow(10.0, -4.0 * u(rng));
|
|
const double I2 = snr * snr * s2;
|
|
const double p = 0.1 + 0.9 * u(rng);
|
|
const double P2 = std::pow((1.0 - p) / p, 2.0) * I2;
|
|
const double e = z(rng);
|
|
pool.push_back({s2, I2, P2, (a * s2 + b * b * I2 + c * c * P2) * e * e});
|
|
}
|
|
const auto fit = FitStillsErrorModel(pool);
|
|
REQUIRE(fit.active);
|
|
const auto many = FitStillsErrorModel(pool, 16); // the same numbers on any number of threads
|
|
CHECK(many.a == fit.a);
|
|
CHECK(many.b2 == fit.b2);
|
|
CHECK(many.c2 == fit.c2);
|
|
CHECK(fit.a == Catch::Approx(a).epsilon(0.05));
|
|
CHECK(std::sqrt(fit.b2) == Catch::Approx(b).epsilon(0.05));
|
|
CHECK(std::sqrt(fit.c2) == Catch::Approx(c).epsilon(0.05));
|
|
}
|
|
|
|
// Without a partiality error in the data the term stays at (near) zero.
|
|
TEST_CASE("StillsErrorModel_NoPartialityErrorNoTerm") {
|
|
std::mt19937 rng(8);
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
std::normal_distribution<double> z(0.0, 1.0);
|
|
std::vector<StillsErrorModelSample> pool;
|
|
for (int i = 0; i < 200000; ++i) {
|
|
const double s2 = std::pow(10.0, 2.0 * u(rng));
|
|
const double snr = 300.0 * std::pow(10.0, -4.0 * u(rng));
|
|
const double I2 = snr * snr * s2;
|
|
const double p = 0.1 + 0.9 * u(rng);
|
|
const double e = z(rng);
|
|
pool.push_back({s2, I2, std::pow((1.0 - p) / p, 2.0) * I2, (1.1 * s2 + 0.0025 * I2) * e * e});
|
|
}
|
|
const auto fit = FitStillsErrorModel(pool);
|
|
REQUIRE(fit.active);
|
|
CHECK(std::sqrt(fit.c2) < 0.01);
|
|
CHECK(std::sqrt(fit.b2) == Catch::Approx(0.05).epsilon(0.05));
|
|
}
|
|
|
|
// 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.b_resolved); // b held at zero is not a b to report 1/b of
|
|
CHECK(fit.b2 == 0.0);
|
|
CHECK(fit.a == Catch::Approx(0.9).epsilon(0.03));
|
|
}
|
|
|
|
// Nearly half the observations repeat their mates to the bit (an input read twice): the bins' median deviations
|
|
// sit near zero, their 1/bd2^2 weights take the fit, and both terms are pulled towards zero. a stops on
|
|
// its lower bound while b^2 goes on to a remainder whose 1/b is an I/sigma in the hundreds or thousands. Beside a
|
|
// clamped a the b^2 is not the solution of any fit, so it is left unresolved and no ISa is printed.
|
|
TEST_CASE("ErrorModel_HalfTheMatesIdenticalLeavesISaUndetermined") {
|
|
std::mt19937 rng(1);
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
std::normal_distribution<double> z(0.0, 1.0);
|
|
std::vector<ErrorModelSample> pool;
|
|
for (int i = 0; i < 200000; ++i) {
|
|
const double s2 = std::pow(10.0, 2.0 * u(rng));
|
|
const double snr = 300.0 * std::pow(10.0, -4.0 * u(rng));
|
|
const double I2 = snr * snr * s2;
|
|
const double e = z(rng);
|
|
double dev2 = (1.3 * s2 + 0.03 * 0.03 * I2) * e * e;
|
|
if (u(rng) < 0.48) dev2 = 0.0;
|
|
pool.push_back({s2, I2, dev2, 2.0f});
|
|
}
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto fit = FitErrorModel(pool, scratch, 4);
|
|
REQUIRE(fit.active);
|
|
CHECK(fit.a == 0.25);
|
|
CHECK(fit.a_at_bound);
|
|
CHECK(1.0 / std::sqrt(fit.b2) > 100.0);
|
|
CHECK(!fit.b_resolved);
|
|
}
|
|
|
|
// The bootstrap over reflections brackets the b the samples were drawn with, and gives the same
|
|
// interval on any number of threads.
|
|
TEST_CASE("ErrorModel_BootstrapIntervalBracketsB") {
|
|
const auto pool = SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.003);
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto fit = FitErrorModel(pool, scratch, 4);
|
|
const auto ci = BootstrapErrorModelB2(pool, fit.b2, 50, pool.size(), 4);
|
|
REQUIRE(ci);
|
|
CHECK(ci->low > 0.0);
|
|
CHECK(ci->low <= 0.03 * 0.03 * 1.1);
|
|
CHECK(ci->high >= 0.03 * 0.03 * 0.9);
|
|
CHECK(ci->high < 2.0 * 0.03 * 0.03);
|
|
const auto one = BootstrapErrorModelB2(pool, fit.b2, 50, pool.size(), 1);
|
|
REQUIRE(one);
|
|
CHECK(one->low == ci->low);
|
|
CHECK(one->high == ci->high);
|
|
// Thinned to a quarter of the pool, the interval comes out about as wide as the full resampling's.
|
|
const auto thin = BootstrapErrorModelB2(pool, fit.b2, 50, pool.size() / 4, 4);
|
|
REQUIRE(thin);
|
|
CHECK(thin->low > 0.0);
|
|
CHECK(thin->high - thin->low == Catch::Approx(ci->high - ci->low).epsilon(0.5));
|
|
}
|
|
|
|
// Data with no systematic error but strong enough for b to be measured: the interval of b^2 reaches
|
|
// zero, so only a lower bound on ISa is measured.
|
|
TEST_CASE("ErrorModel_BootstrapIntervalReachesZeroWithoutB") {
|
|
const auto pool = SyntheticErrorModelSamples(1.3, 0.0, 300.0, 0.0);
|
|
std::vector<ErrorModelBinned> scratch;
|
|
const auto fit = FitErrorModel(pool, scratch, 4);
|
|
REQUIRE(fit.b_measured);
|
|
const auto ci = BootstrapErrorModelB2(pool, fit.b2, 50, pool.size(), 4);
|
|
REQUIRE(ci);
|
|
CHECK(ci->low == 0.0);
|
|
CHECK(ci->high > 0.0);
|
|
}
|
|
|
|
#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_split == parallel_stats.overall.r_split);
|
|
// R_merge / R_meas / R_pim are rotation statistics: the stills merge does not form them.
|
|
CHECK(std::isnan(serial_stats.overall.r_meas));
|
|
CHECK(std::isnan(serial_stats.overall.r_merge));
|
|
CHECK(std::isnan(serial_stats.overall.r_pim));
|
|
}
|
|
|
|
// With the error model fitted, the merge takes every observation, however small its partiality, but reports a
|
|
// reflection only if one of its observations reaches the existence partiality (0.2).
|
|
TEST_CASE("MergeOnTheFly_EveryObservationMergedReflectionReportedOnlyIfOneReachesTheExistencePartiality") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(TestSpaceGroup());
|
|
ScalingSettings s = x.GetScalingSettings();
|
|
x.ImportScalingSettings(s.MinPartiality(0.5)); // the stills default
|
|
|
|
std::mt19937 rng(7);
|
|
std::uniform_int_distribution<int> index(-6, 6);
|
|
std::uniform_real_distribution<float> unit(0.0f, 1.0f);
|
|
auto observation = [](int h, int k, int l, float I, float sigma, float p) {
|
|
Reflection r{};
|
|
r.h = h;
|
|
r.k = k;
|
|
r.l = l;
|
|
r.d = 20.0f / std::sqrt(static_cast<float>(h * h + k * k + l * l));
|
|
r.I = I;
|
|
r.sigma = sigma;
|
|
r.var_bkg = 0.5f * sigma * sigma;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = p;
|
|
r.image_scale_corr = 1.0f;
|
|
return r;
|
|
};
|
|
// Enough ordinary observations for the error model.
|
|
std::vector<IntegrationOutcome> outcomes(5000);
|
|
for (auto &o: outcomes)
|
|
for (int j = 0; j < 20; ++j) {
|
|
const int h = index(rng), k = index(rng), l = index(rng);
|
|
if (h != 0 || k != 0 || l != 0)
|
|
o.reflections.push_back(observation(h, k, l, 100.0f * unit(rng) - 5.0f, 1.0f + 5.0f * unit(rng),
|
|
0.5f + 0.5f * unit(rng)));
|
|
}
|
|
// (8,1,1): seen only in the tails. (9,1,1): one observation at p = 0.3 and three at p = 0.1.
|
|
// (9,2,1): the same p = 0.3 observation alone.
|
|
for (int i = 0; i < 4; ++i)
|
|
outcomes[i].reflections.push_back(observation(8, 1, 1, 50.0f, 2.0f, 0.1f));
|
|
outcomes[10].reflections.push_back(observation(9, 1, 1, 50.0f, 2.0f, 0.3f));
|
|
for (int i = 11; i < 14; ++i)
|
|
outcomes[i].reflections.push_back(observation(9, 1, 1, 50.0f, 2.0f, 0.1f));
|
|
outcomes[20].reflections.push_back(observation(9, 2, 1, 50.0f, 2.0f, 0.3f));
|
|
|
|
MergeOnTheFly merge(x);
|
|
merge.RefineErrorModel(outcomes);
|
|
REQUIRE(merge.ErrorModelActive());
|
|
merge.AddImages(outcomes);
|
|
const auto merged = merge.ExportReflections();
|
|
auto find = [&](int h, int k, int l) -> const MergedReflection * {
|
|
for (const auto &m: merged)
|
|
if (m.h == h && m.k == k && m.l == l)
|
|
return &m;
|
|
return nullptr;
|
|
};
|
|
CHECK(find(8, 1, 1) == nullptr);
|
|
const auto *with_tails = find(9, 1, 1);
|
|
const auto *alone = find(9, 2, 1);
|
|
REQUIRE(with_tails != nullptr);
|
|
REQUIRE(alone != nullptr);
|
|
CHECK(with_tails->sigma < alone->sigma); // the p = 0.1 observations are in its merge
|
|
|
|
// None of the tail observations is flagged as below the merge's floor.
|
|
const auto verdicts = merge.Verdicts(outcomes);
|
|
CHECK((verdicts[0].back().flags & OBS_LOW_PARTIALITY) == 0);
|
|
CHECK((verdicts[11].back().flags & OBS_LOW_PARTIALITY) == 0);
|
|
}
|
|
|
|
// The existence rule holds within each half-set too: a half that has only tail observations of a reflection has no
|
|
// mean of it. And observations far out on the rocking curve (p = 1e-30, a correction no float variance holds)
|
|
// leave every merged value and statistic finite.
|
|
TEST_CASE("MergeOnTheFly_HalfSetWithOnlyTailObservationsHasNoMean") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(TestSpaceGroup());
|
|
ScalingSettings s = x.GetScalingSettings();
|
|
x.ImportScalingSettings(s.MinPartiality(0.5));
|
|
|
|
std::mt19937 rng(7);
|
|
std::uniform_int_distribution<int> index(-6, 6);
|
|
std::uniform_real_distribution<float> unit(0.0f, 1.0f);
|
|
auto observation = [](int h, int k, int l, float I, float sigma, float p, float corr) {
|
|
Reflection r{};
|
|
r.h = h;
|
|
r.k = k;
|
|
r.l = l;
|
|
r.d = 20.0f / std::sqrt(static_cast<float>(h * h + k * k + l * l));
|
|
r.I = I;
|
|
r.sigma = sigma;
|
|
r.var_bkg = 0.5f * sigma * sigma;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = p;
|
|
r.image_scale_corr = corr;
|
|
return r;
|
|
};
|
|
std::vector<IntegrationOutcome> outcomes(5000);
|
|
for (auto &o: outcomes)
|
|
for (int j = 0; j < 20; ++j) {
|
|
const int h = index(rng), k = index(rng), l = index(rng);
|
|
if (h != 0 || k != 0 || l != 0)
|
|
o.reflections.push_back(observation(h, k, l, 100.0f * unit(rng) - 5.0f, 1.0f + 5.0f * unit(rng),
|
|
0.5f + 0.5f * unit(rng), 1.0f));
|
|
}
|
|
// Far-tail observations of the ordinary reflections, with the corrections such a partiality brings.
|
|
for (size_t i = 0; i < outcomes.size(); i += 5)
|
|
outcomes[i].reflections.push_back(observation(1, 2, 3, 1e4f, 30.0f, 1e-30f, 1e36f));
|
|
// (9,1,1): one observation at p = 0.3 on image 10, then tail observations on twenty more images - the
|
|
// other half-set, and the other frame-whole half-set, have only tail observations.
|
|
outcomes[10].reflections.push_back(observation(9, 1, 1, 50.0f, 2.0f, 0.3f, 1.0f));
|
|
for (int i = 11; i < 31; ++i)
|
|
outcomes[i].reflections.push_back(observation(9, 1, 1, 0.1f, 2.0f, 0.01f, 100.0f));
|
|
|
|
MergeOnTheFly merge(x);
|
|
merge.RefineErrorModel(outcomes);
|
|
REQUIRE(merge.ErrorModelActive());
|
|
merge.AddImages(outcomes);
|
|
CHECK(merge.FarTailCount() == 1000);
|
|
const auto merged = merge.ExportReflections();
|
|
|
|
const MergedReflection *m = nullptr;
|
|
for (const auto &r: merged) {
|
|
if (r.h == 9 && r.k == 1 && r.l == 1)
|
|
m = &r;
|
|
CHECK(std::isfinite(r.I));
|
|
CHECK(std::isfinite(r.sigma));
|
|
CHECK(r.sigma > 0.0f);
|
|
}
|
|
REQUIRE(m != nullptr);
|
|
CHECK(!std::isfinite(m->I_half[0]));
|
|
CHECK(!std::isfinite(m->I_half[1]));
|
|
CHECK(std::isfinite(m->I_split[0]) != std::isfinite(m->I_split[1]));
|
|
|
|
const auto stats = merge.MergeStats(merged, outcomes);
|
|
CHECK(std::isfinite(stats.overall.cc_half));
|
|
CHECK(std::isfinite(stats.overall.cc_half_split));
|
|
CHECK(std::isfinite(stats.overall.r_split));
|
|
}
|
|
|
|
// 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.dist_ewald = std::fabs((q + S0).Length() - 1.0f / WAVELENGTH);
|
|
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.41 (0.54 with the per-crystal window width): 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 reciprocal lattice scaled by (1 + eps) - a wavelength or cell a little off, which the spot positions
|
|
// cannot see - centres every crystal's rocking curve at eps lambda d*^2 / 2 instead of zero. The windows are
|
|
// cut symmetrically about the stored lattice's sphere, as prediction cuts them; the run's radial scale has
|
|
// to come back from where the recorded intensity sits within them.
|
|
TEST_CASE("StillsPartialityRefine_MeasureRadialScale") {
|
|
constexpr float CELL = 50.0f, WAVELENGTH = 1.0f, D_MIN = 2.0f;
|
|
constexpr double EPS = 4e-3, GAMMA0 = 2e-4, WINDOW = 1e-3;
|
|
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(13);
|
|
std::normal_distribution<double> gauss(0.0, 1.0);
|
|
std::uniform_real_distribution<double> unit(0.0, 1.0);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::vector<IntegrationOutcome> outcomes(100);
|
|
const int hmax = static_cast<int>(CELL / D_MIN);
|
|
for (auto &o: outcomes) {
|
|
const auto axes = TurnedAxes(UnitCell(CELL, CELL, CELL, 90, 90, 90));
|
|
// a different orientation per crystal: turn the turned axes about a random axis
|
|
const double ax = gauss(rng), ay = gauss(rng), az = gauss(rng), n = std::sqrt(ax * ax + ay * ay + az * az);
|
|
const double angle = 2.0 * M_PI * unit(rng), c = std::cos(angle), s = std::sin(angle);
|
|
auto turn = [&](const Coord &v) {
|
|
const double k[3] = {ax / n, ay / n, az / n};
|
|
const double dot = k[0] * v.x + k[1] * v.y + k[2] * v.z;
|
|
const double cr[3] = {k[1] * v.z - k[2] * v.y, k[2] * v.x - k[0] * v.z, k[0] * v.y - k[1] * v.x};
|
|
return Coord(static_cast<float>(v.x * c + cr[0] * s + k[0] * dot * (1 - c)),
|
|
static_cast<float>(v.y * c + cr[1] * s + k[1] * dot * (1 - c)),
|
|
static_cast<float>(v.z * c + cr[2] * s + k[2] * dot * (1 - c)));
|
|
};
|
|
o.geom = geom;
|
|
o.latt = CrystalLattice(turn(axes[0]), turn(axes[1]), turn(axes[2]));
|
|
o.image_scale_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;
|
|
const double de_stored = (q + S0).Length() - 1.0 / WAVELENGTH;
|
|
if (std::fabs(de_stored) > WINDOW)
|
|
continue;
|
|
const double de_true = (q * static_cast<float>(1.0 + EPS) + S0).Length() - 1.0 / WAVELENGTH;
|
|
Reflection r{};
|
|
r.h = h; r.k = k; r.l = l;
|
|
r.d = 1.0f / dstar;
|
|
r.dist_ewald = static_cast<float>(std::fabs(de_stored));
|
|
r.I = static_cast<float>(*o.image_scale_g * 1000.0 * wilson(rng)
|
|
* std::exp(-0.5 * de_true * de_true / (GAMMA0 * GAMMA0)));
|
|
r.sigma = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 1.0f;
|
|
o.reflections.push_back(r);
|
|
}
|
|
}
|
|
const StillsPartialityRefine refine(x);
|
|
double scale = 1.0;
|
|
for (int i = 0; i < 3; ++i)
|
|
scale = refine.MeasureRadialScale(outcomes, std::nullopt, scale);
|
|
CHECK(scale - 1.0 == Catch::Approx(EPS).epsilon(0.15));
|
|
|
|
// and a lattice that is right stays put
|
|
for (auto &o: outcomes)
|
|
for (auto &r: o.reflections)
|
|
r.I = static_cast<float>(*o.image_scale_g * 1000.0);
|
|
CHECK(std::fabs(refine.MeasureRadialScale(outcomes, std::nullopt, 1.0) - 1.0) < 5e-4);
|
|
}
|
|
|
|
|
|
// 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);
|
|
}
|
|
}
|
|
|
|
// A cell built on face diagonals (a+b, a-b, c) indexes every spot of the true cell, all at h+k even: the
|
|
// sub-lattice that drops the h+k odd nodes is taken, and its cell is the true one. Where those nodes hold
|
|
// spots like the others - a real superstructure - or there are too few spots to tell, the cell is kept.
|
|
TEST_CASE("TestSubLattices_FaceDiagonalSupercellComesDownToTheTrueCell") {
|
|
const UnitCell primitive{70, 71, 73, 90, 90, 90}; // synthetic
|
|
const CrystalLattice p(primitive);
|
|
const CrystalLattice super(p.Vec0() + p.Vec1(), p.Vec0() - p.Vec1(), p.Vec2());
|
|
const UnitCell super_cell = super.GetUnitCell();
|
|
|
|
const auto frame = [](bool superstructure, int n) {
|
|
FrameSpotIndices f;
|
|
for (int i = 0; i < n; i++) {
|
|
const int h = i % 7 - 3, k = (i / 7) % 7 - 3, l = i / 49 - 2;
|
|
if (superstructure || (h + k) % 2 == 0)
|
|
f.hkl.push_back({h, k, l});
|
|
}
|
|
f.off_lattice = 40;
|
|
return f;
|
|
};
|
|
const FrameSpotIndices empty_between = frame(false, 245), real_between = frame(true, 245), few = frame(false, 6);
|
|
|
|
const auto t = TestSubLattices({&empty_between, &empty_between}, 0.1f);
|
|
CHECK(t.adopt);
|
|
CHECK(t.index == 2);
|
|
CHECK(t.dropped == 0);
|
|
CHECK(t.translation == std::array<int, 3>{1, 1, 0});
|
|
const UnitCell sub = SubLatticeCell(super_cell, t.translation, t.index);
|
|
CHECK(std::abs(std::abs(CrystalLattice(sub).CalcVolume()) - std::abs(p.CalcVolume())) < 1.0);
|
|
CHECK(sub.is_close(primitive, 0.001, 0.01));
|
|
|
|
CHECK_FALSE(TestSubLattices({&real_between, &real_between}, 0.1f).adopt);
|
|
CHECK_FALSE(TestSubLattices({&few}, 0.1f).adopt);
|
|
}
|
|
|
|
// R_split's half-sets keep every frame whole: each frame measures the same reflections at its own
|
|
// intensity, so every reflection's two half-set means are the same two numbers. R_split is then
|
|
// (1/sqrt 2) |I_1 - I_2| / ((I_1 + I_2) / 2) in every shell, and the frame-split CC1/2 is undefined
|
|
// (the half-set means do not vary) - the per-reflection balanced split, which the processing decides on,
|
|
// is left as it was.
|
|
TEST_CASE("MergeOnTheFly_RsplitKeepsEveryFrameWhole") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(TestSpaceGroup());
|
|
|
|
std::vector<IntegrationOutcome> outcomes(64);
|
|
for (size_t f = 0; f < outcomes.size(); ++f) {
|
|
for (int h = 1; h <= 4; ++h)
|
|
for (int k = 1; k <= 4; ++k) {
|
|
Reflection r{};
|
|
r.h = h;
|
|
r.k = k;
|
|
r.l = 1;
|
|
r.d = 20.0f / std::sqrt(static_cast<float>(h * h + k * k + 1));
|
|
r.I = 100.0f + static_cast<float>(f);
|
|
r.sigma = 10.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 1.0f;
|
|
r.image_scale_corr = 1.0f;
|
|
outcomes[f].reflections.push_back(r);
|
|
}
|
|
}
|
|
|
|
MergeOnTheFly merge(x);
|
|
merge.AddImages(outcomes);
|
|
const auto merged = merge.ExportReflections();
|
|
REQUIRE(merged.size() == 16);
|
|
|
|
const float a = merged[0].I_split[0], b = merged[0].I_split[1];
|
|
REQUIRE(std::isfinite(a));
|
|
REQUIRE(std::isfinite(b));
|
|
CHECK(a != b);
|
|
for (const auto &m: merged) {
|
|
CHECK(m.I_split[0] == a);
|
|
CHECK(m.I_split[1] == b);
|
|
}
|
|
|
|
const auto stats = merge.MergeStats(merged, outcomes);
|
|
const double expected = std::fabs(a - b) / (0.5 * (a + b)) / std::sqrt(2.0);
|
|
CHECK(stats.overall.n_split_pairs == 16);
|
|
CHECK(stats.overall.r_split == Catch::Approx(expected).epsilon(1e-5));
|
|
for (const auto &sh: stats.shells)
|
|
if (sh.n_split_pairs > 0)
|
|
CHECK(sh.r_split == Catch::Approx(expected).epsilon(1e-5));
|
|
CHECK_FALSE(SerialShellTable(stats).empty());
|
|
}
|
|
|
|
// The run's rocking width is measured on the crystals, not on the prediction window: reflections spread
|
|
// evenly across a window four times wider than the curve, each recording I_full * p(de), give back the
|
|
// curve's gamma0 and angular width, while the plain moment of de would return the window's.
|
|
TEST_CASE("StillsPartialityRefine_MeasureRunWidth") {
|
|
constexpr double GAMMA0 = 1.0e-4, GAMMA_E = 0.035 * M_PI / 180.0;
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
|
|
std::mt19937 rng(5);
|
|
std::uniform_real_distribution<double> unit(0.0, 1.0);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::normal_distribution<double> gauss(0.0, 1.0);
|
|
std::vector<IntegrationOutcome> outcomes(200);
|
|
for (auto &o: outcomes) {
|
|
o.image_scale_g = 0.5 + unit(rng);
|
|
for (int j = 0; j < 200; ++j) {
|
|
const double d = 1.5 + 8.0 * unit(rng);
|
|
const double dstar = 1.0 / d;
|
|
const double s = std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E * dstar * dstar);
|
|
const double window = 4.0 * std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E / (1.5 * 1.5));
|
|
const double de = window * (2.0 * unit(rng) - 1.0);
|
|
Reflection r{};
|
|
r.h = j; r.k = 1; r.l = 1;
|
|
r.d = static_cast<float>(d);
|
|
r.dist_ewald = static_cast<float>(std::fabs(de));
|
|
r.I = static_cast<float>(*o.image_scale_g * 1000.0 * wilson(rng) * std::exp(-0.5 * de * de / (s * s))
|
|
+ 5.0 * gauss(rng));
|
|
r.sigma = 5.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 1.0f;
|
|
r.image_scale_corr = static_cast<float>(1.0 / *o.image_scale_g);
|
|
o.reflections.push_back(r);
|
|
}
|
|
}
|
|
const auto width = StillsPartialityRefine(x).MeasureRunWidth(outcomes, std::nullopt);
|
|
REQUIRE(width.has_value());
|
|
CHECK(width->gamma0 == Catch::Approx(GAMMA0).epsilon(0.3));
|
|
CHECK(width->gamma_e == Catch::Approx(GAMMA_E).epsilon(0.1));
|
|
CHECK(width->bandwidth_sigma < 1e-3);
|
|
}
|
|
|
|
// Frames predicted through a much wider window whose intensities do not follow the curve (a wrong lattice
|
|
// on the frame) must not decide the measured width: pooled with the others they hold most of the de^2, and
|
|
// the run's width then read as their window. Measured per group of similar window, they are a group whose
|
|
// curve looks wider than its window, which is left out.
|
|
TEST_CASE("StillsPartialityRefine_MeasureRunWidth_WideWindowFramesIgnored") {
|
|
constexpr double GAMMA0 = 1.0e-4, GAMMA_E = 0.035 * M_PI / 180.0;
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
|
|
std::mt19937 rng(7);
|
|
std::uniform_real_distribution<double> unit(0.0, 1.0);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::normal_distribution<double> gauss(0.0, 1.0);
|
|
std::vector<IntegrationOutcome> outcomes(400);
|
|
for (size_t f = 0; f < outcomes.size(); ++f) {
|
|
auto &o = outcomes[f];
|
|
const bool wide = f % 4 == 3;
|
|
o.image_scale_g = 0.5 + unit(rng);
|
|
for (int j = 0; j < 200; ++j) {
|
|
const double d = 1.5 + 8.0 * unit(rng);
|
|
const double dstar = 1.0 / d;
|
|
const double s = std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E * dstar * dstar);
|
|
const double window = (wide ? 20.0 : 4.0) * std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E / (1.5 * 1.5));
|
|
const double de = window * (2.0 * unit(rng) - 1.0);
|
|
const double p = wide ? 1.0 : std::exp(-0.5 * de * de / (s * s));
|
|
Reflection r{};
|
|
r.h = j; r.k = 1; r.l = 1;
|
|
r.d = static_cast<float>(d);
|
|
r.dist_ewald = static_cast<float>(std::fabs(de));
|
|
r.I = static_cast<float>(*o.image_scale_g * 1000.0 * wilson(rng) * p + 5.0 * gauss(rng));
|
|
r.sigma = 5.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 1.0f;
|
|
r.image_scale_corr = static_cast<float>(1.0 / *o.image_scale_g);
|
|
o.reflections.push_back(r);
|
|
}
|
|
}
|
|
const auto width = StillsPartialityRefine(x).MeasureRunWidth(outcomes, std::nullopt);
|
|
REQUIRE(width.has_value());
|
|
CHECK(width->gamma0 == Catch::Approx(GAMMA0).epsilon(0.3));
|
|
CHECK(width->gamma_e == Catch::Approx(GAMMA_E).epsilon(0.15));
|
|
}
|
|
|
|
// Frames predicted through the run's width (2 sigma(d*) of its formula) have a window at every reflection's
|
|
// own resolution; measured with that window, the table gives the width back at each shell.
|
|
TEST_CASE("StillsPartialityRefine_MeasureRunWidth_WindowAtReflectionResolution") {
|
|
constexpr double GAMMA0 = 1.0e-4, GAMMA_E = 0.035 * M_PI / 180.0;
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
const StillsPartialityRefine::RunWidth window{.gamma0 = GAMMA0, .gamma_e = GAMMA_E};
|
|
auto sigma = [&](double dstar) { return std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E * dstar * dstar); };
|
|
|
|
std::mt19937 rng(11);
|
|
std::uniform_real_distribution<double> unit(0.0, 1.0);
|
|
std::exponential_distribution<double> wilson(1.0);
|
|
std::normal_distribution<double> gauss(0.0, 1.0);
|
|
std::vector<IntegrationOutcome> outcomes(400);
|
|
for (auto &o: outcomes) {
|
|
o.image_scale_g = 0.5 + unit(rng);
|
|
for (int j = 0; j < 200; ++j) {
|
|
const double d = 1.5 + 8.0 * unit(rng);
|
|
const double s = sigma(1.0 / d);
|
|
const double de = 2.0 * s * (2.0 * unit(rng) - 1.0);
|
|
Reflection r{};
|
|
r.h = j; r.k = 1; r.l = 1;
|
|
r.d = static_cast<float>(d);
|
|
r.dist_ewald = static_cast<float>(std::fabs(de));
|
|
r.I = static_cast<float>(*o.image_scale_g * 1000.0 * wilson(rng) * std::exp(-0.5 * de * de / (s * s))
|
|
+ 5.0 * gauss(rng));
|
|
r.sigma = 5.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 1.0f;
|
|
r.image_scale_corr = static_cast<float>(1.0 / *o.image_scale_g);
|
|
o.reflections.push_back(r);
|
|
}
|
|
}
|
|
const auto width = StillsPartialityRefine(x).MeasureRunWidth(outcomes, std::nullopt, 1.0, &window);
|
|
REQUIRE(width.has_value());
|
|
REQUIRE(width->table_s.size() >= 3);
|
|
for (size_t i = 0; i < width->table_s.size(); ++i) {
|
|
INFO("d* = " << width->table_dstar[i]);
|
|
CHECK(width->table_s[i] == Catch::Approx(sigma(width->table_dstar[i])).epsilon(0.15));
|
|
}
|
|
}
|
|
|
|
// Outlier rejection on the stills merge: one observation at ten times its reflection's intensity among
|
|
// nineteen consistent ones is dropped, and nothing else is.
|
|
TEST_CASE("MergeOnTheFly_RejectsAnOutlierInTheBand") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
ScalingSettings s = x.GetScalingSettings();
|
|
s.OutlierRejectNsigma(6.0);
|
|
x.ImportScalingSettings(s);
|
|
|
|
std::mt19937 rng(3);
|
|
std::normal_distribution<float> gauss(0.0f, 1.0f);
|
|
std::vector<IntegrationOutcome> outcomes(20);
|
|
for (size_t i = 0; i < outcomes.size(); ++i) {
|
|
for (int h = 1; h <= 50; ++h) {
|
|
Reflection r{};
|
|
r.h = h; r.k = 2; r.l = 3;
|
|
r.d = 2.0f;
|
|
const float I_true = 100.0f * static_cast<float>(h);
|
|
r.I = I_true * (1.0f + 0.1f * gauss(rng)) + 3.0f * gauss(rng);
|
|
if (i == 7 && h == 25)
|
|
r.I = 10.0f * I_true;
|
|
r.sigma = 3.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 1.0f;
|
|
r.image_scale_corr = 1.0f;
|
|
outcomes[i].reflections.push_back(r);
|
|
}
|
|
}
|
|
MergeOnTheFly merge(x);
|
|
merge.RefineErrorModel(outcomes);
|
|
REQUIRE(merge.ErrorModelActive());
|
|
merge.AddImages(outcomes);
|
|
CHECK(merge.RejectedCount() == 1);
|
|
}
|
|
|
|
// Stills: the observations of a reflection spread from near zero to its full intensity (partials), so
|
|
// the spread of ln(I / median) is wide. One powder grain in a weak reflection's box, at 200 times its
|
|
// intensity, must still be rejected - the band is the error model's, n sigma about the median - while
|
|
// none of the partials is.
|
|
TEST_CASE("MergeOnTheFly_StillsRejectsAGrainAmongPartials") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
ScalingSettings s = x.GetScalingSettings();
|
|
s.OutlierRejectNsigma(6.0);
|
|
x.ImportScalingSettings(s);
|
|
REQUIRE(!x.IsRotationIndexing());
|
|
|
|
std::mt19937 rng(5);
|
|
std::normal_distribution<float> gauss(0.0f, 1.0f);
|
|
std::uniform_real_distribution<float> partial(0.1f, 1.0f);
|
|
std::vector<IntegrationOutcome> outcomes(40);
|
|
for (size_t i = 0; i < outcomes.size(); ++i) {
|
|
for (int h = 1; h <= 50; ++h) {
|
|
Reflection r{};
|
|
r.h = h; r.k = 2; r.l = 3;
|
|
r.d = 2.0f;
|
|
const float I_true = 100.0f * static_cast<float>(h);
|
|
r.I = I_true * partial(rng) + 3.0f * gauss(rng);
|
|
if (i == 7 && h == 1)
|
|
r.I = 200.0f * I_true;
|
|
r.sigma = 3.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.partiality = 1.0f;
|
|
r.image_scale_corr = 1.0f;
|
|
outcomes[i].reflections.push_back(r);
|
|
}
|
|
}
|
|
MergeOnTheFly merge(x);
|
|
merge.RefineErrorModel(outcomes);
|
|
REQUIRE(merge.ErrorModelActive());
|
|
merge.AddImages(outcomes);
|
|
CHECK(merge.RejectedCount() == 1);
|
|
}
|
|
|
|
// The capture term widens the sigma of a partial by (c (1 - p) <I>): a partial merges at a lower weight
|
|
// than a full of the same intensity, and not at all differently with c = 0.
|
|
TEST_CASE("MergeOnTheFly_CaptureTermDownweightsPartials") {
|
|
auto merged_I = [](double c) {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
ScalingSettings s = x.GetScalingSettings();
|
|
s.CaptureUncertaintyCoeff(c);
|
|
x.ImportScalingSettings(s);
|
|
std::mt19937 rng(9);
|
|
std::normal_distribution<float> gauss(0.0f, 1.0f);
|
|
std::vector<IntegrationOutcome> outcomes(40);
|
|
for (size_t i = 0; i < outcomes.size(); ++i)
|
|
for (int h = 1; h <= 40; ++h) {
|
|
Reflection r{};
|
|
r.h = h; r.k = 1; r.l = 1;
|
|
r.d = 2.0f;
|
|
// Every other image records the reflection as a 30 % partial that reads 20 % high.
|
|
const bool partial = i % 2 == 1;
|
|
r.partiality = partial ? 0.3f : 1.0f;
|
|
r.image_scale_corr = 1.0f / r.partiality;
|
|
r.I = r.partiality * 100.0f * static_cast<float>(h) * (partial ? 1.2f : 1.0f) * (1.0f + 0.05f * gauss(rng));
|
|
r.sigma = 2.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
outcomes[i].reflections.push_back(r);
|
|
}
|
|
MergeOnTheFly merge(x);
|
|
merge.RefineErrorModel(outcomes);
|
|
merge.AddImages(outcomes);
|
|
const auto m = merge.ExportReflections();
|
|
return m.at(10).I / (100.0f * static_cast<float>(m.at(10).h));
|
|
};
|
|
const float without = merged_I(0.0), with = merged_I(1.0);
|
|
CHECK(with < without);
|
|
CHECK(with < 1.08f);
|
|
}
|
|
|
|
// The setter checks the value it is given, not the one it replaces.
|
|
TEST_CASE("ScalingSettings_MinPartialityRange") {
|
|
ScalingSettings s;
|
|
CHECK_THROWS(s.MinPartiality(1.5));
|
|
CHECK_THROWS(s.MinPartiality(-0.1));
|
|
CHECK_NOTHROW(s.MinPartiality(0.5));
|
|
CHECK(s.GetMinPartiality() == 0.5);
|
|
}
|
|
|
|
// An image with too few reflections to fit its scale on is left out of the merge, not merged at the scale
|
|
// of 1 it was predicted with.
|
|
TEST_CASE("ScaleOnTheFly_TooFewReflectionsLeavesTheImageOut") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
std::vector<IntegrationOutcome> outcomes(3);
|
|
for (size_t i = 0; i < outcomes.size(); ++i) {
|
|
const int n = i == 2 ? 10 : 50;
|
|
for (int h = 1; h <= n; ++h) {
|
|
Reflection r{};
|
|
r.h = h; r.k = 1; r.l = 1;
|
|
r.d = 3.0f;
|
|
r.I = 100.0f * static_cast<float>(h) * (i == 1 ? 2.0f : 1.0f);
|
|
r.sigma = 1.0f;
|
|
r.var_bkg = 0.5f;
|
|
r.partiality = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.image_scale_corr = 1.0f;
|
|
outcomes[i].reflections.push_back(r);
|
|
}
|
|
}
|
|
ScaleOnTheFly(x, MergeAll(x, outcomes, 1, std::nullopt)).Scale(outcomes, 1);
|
|
CHECK(outcomes[0].image_scale_g.has_value());
|
|
CHECK(std::isfinite(outcomes[1].reflections[0].image_scale_corr));
|
|
CHECK(!outcomes[2].image_scale_g.has_value());
|
|
for (const auto &r: outcomes[2].reflections)
|
|
CHECK(std::isnan(r.image_scale_corr));
|
|
}
|
|
|
|
// A frame whose background region reads below zero (a module not counting photons) integrates reflections with
|
|
// a negative background variance, a sigma at its floor of one photon and an intensity of thousands. Merged,
|
|
// such an observation takes its reflection over, and every frame that sees the reflection - and finds
|
|
// nothing there - is then fitted towards a scale of 0. Left out, the frames keep their own scales.
|
|
TEST_CASE("ScaleOnTheFly_NegativeBackgroundVarianceDoesNotCollapseTheScales") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
constexpr int frames = 20, n = 60;
|
|
std::vector<IntegrationOutcome> outcomes(frames);
|
|
for (int i = 0; i < frames; ++i) {
|
|
const float G = 1.0f + 0.1f * static_cast<float>(i % 5);
|
|
for (int h = 1; h <= n; ++h) {
|
|
Reflection r{};
|
|
r.h = h; r.k = 1; r.l = 1;
|
|
r.d = 3.0f;
|
|
// the first ten are beyond the crystal's diffraction: nothing there
|
|
r.I = h <= 10 ? 0.0f : G * 100.0f * static_cast<float>(1 + h % 7);
|
|
r.var_bkg = 50.0f;
|
|
r.sigma = std::sqrt(r.I + r.var_bkg);
|
|
r.partiality = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.image_scale_corr = 1.0f;
|
|
if (i == 0 && h <= 10) { // the broken module of frame 0
|
|
r.I = 1e5f;
|
|
r.var_bkg = -1e4f;
|
|
r.sigma = 1.0f;
|
|
}
|
|
outcomes[i].reflections.push_back(r);
|
|
}
|
|
}
|
|
ScaleOnTheFly(x, MergeAll(x, outcomes, 1, std::nullopt)).Scale(outcomes, 1);
|
|
for (int i = 1; i < frames; ++i) {
|
|
REQUIRE(outcomes[i].image_scale_g.has_value());
|
|
CHECK(*outcomes[i].image_scale_g / *outcomes[1].image_scale_g
|
|
== Catch::Approx((1.0 + 0.1 * (i % 5)) / 1.1).epsilon(0.01));
|
|
}
|
|
}
|
|
|
|
// Stills with the Friedel hands kept apart report CC_anom (between frame-whole half-sets) and SigAno: a
|
|
// consistent anomalous difference far above the noise correlates close to 1.
|
|
TEST_CASE("MergeOnTheFly_StillsAnomalousStatistics") {
|
|
DiffractionExperiment x;
|
|
x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
|
|
ScalingSettings s = x.GetScalingSettings();
|
|
s.MergeFriedel(false);
|
|
x.ImportScalingSettings(s);
|
|
std::mt19937 rng(11);
|
|
std::normal_distribution<float> gauss(0.0f, 1.0f);
|
|
std::vector<float> dano(61);
|
|
for (auto &v: dano) v = 0.2f * gauss(rng);
|
|
std::vector<IntegrationOutcome> outcomes(40);
|
|
for (size_t i = 0; i < outcomes.size(); ++i)
|
|
for (int h = 1; h <= 60; ++h)
|
|
for (int sign: {1, -1}) {
|
|
Reflection r{};
|
|
r.h = sign * h; r.k = sign * 2; r.l = sign * 3;
|
|
r.d = 30.0f / std::sqrt(static_cast<float>(h * h + 13));
|
|
r.I = 1000.0f * (1.0f + 0.5f * sign * dano[h]) + 5.0f * gauss(rng);
|
|
r.sigma = 5.0f;
|
|
r.partiality = 1.0f;
|
|
r.prescaling_corr = 1.0f;
|
|
r.image_scale_corr = 1.0f;
|
|
outcomes[i].reflections.push_back(r);
|
|
}
|
|
MergeOnTheFly merge(x);
|
|
merge.AddImages(outcomes);
|
|
const auto merged = merge.ExportReflections();
|
|
const auto stats = merge.MergeStats(merged, outcomes);
|
|
CHECK(stats.overall.cc_anom > 0.9);
|
|
CHECK(stats.overall.abs_diff_over_sigma_anomalous > 5.0);
|
|
}
|