mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity
Files
02ea56b2d320d396103bc2d52182451527fef94a
Jungfraujoch/tests/MergeScaleTest.cpp
T
leonarski_fandClaude Opus 5.5 3c4d730f89 RESULT-CHANGING (owner decision): stills tilt post-refinement solved directly instead of by Ceres 
Each crystal's (psi_x, psi_y) tilt was refined by a Ceres problem rebuilt for every inner
iteration: one auto-differentiated residual block and one Cauchy loss per reflection. With the
scaling reference cleaned up it was 70% of the stills scale/merge. The same objective - Cauchy
on the sigma-weighted intensity residuals plus the Gaussian prior on the tilt, the tilt bounded -
is now minimised by a small Levenberg-Marquardt on the rho'-reweighted 2x2 normal equations,
the tilt applied to all reflections through one rotation matrix and the Jacobian taken by
central differences.

It reaches the same answer to optimizer precision: on the new synthetic test
(StillsPartialityRefine_RecoversTheTilt) partiality error 0.0559 -> 0.03019 with Ceres and
0.03016 with this, mean tilt 0.02065 vs 0.02067 deg; on two 25k-frame XFEL runs (--mode scale)
CC1/2, I/sigma, R_pim and the CC against a deposited lysozyme model agree to 0.0005, mean
|dpsi| 0.023 deg on both. Scale+merge 6.9 -> 3.0 s (--mode scale, loaded box).

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
2026-10-14 11:59:09 +02:00

749 lines
36 KiB
C++
Raw Blame History

// 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 "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], &parallel_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
}
Reference in New Issue View Git Blame Copy Permalink