mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases 73 Wiki Activity
Files
0c3e1462ef351b678ad0fb0d060e6b032ed815d1
Jungfraujoch/tests/MergeScaleTest.cpp
T
leonarski_fandClaude Opus 5 0c3e1462ef merge: the completeness denominator is the declared range, not the surviving one 
CalcPossibleReflections was handed d_min/d_max derived from the reflections that
came out of the merge, so a loss at either extreme took the numerator and the
denominator with it. At the high end that is right: d_min is the finest d reached
anywhere and the denominator is the full sphere down to it, so anisotropic loss
shows. At the low end it was a tautology - d_max was the coarsest reflection that
happened to survive, so anything the beam-stop shadow mask (on by default), a
detector mask or the low-resolution limit itself removed left the denominator
along with the data and could not be reported as missing.

Both statistics paths now bin, and count, between the DECLARED low-resolution
limit and the finest d reached: MergeOnTheFly::MergeStats (stills) and
RotationScaleMerge::MergeAndStats (rotation). The grid and the denominator keep
sharing their bounds, so no possible reflection falls outside a shell. An
undeclared low limit is the whole sphere - 1/d^2 down to 0 - spelled as an
infinite d_max, which ResolutionShells already handles and which gemmi's
for_all_reflections special-cases; the change therefore reads correctly whether
or not the 50 A default stays. The innermost shell keeps a finite d_max label,
falling back to the coarsest reflection measured when the bound is infinite.

This makes the shell boundaries the ones the integration document already claims:
XDS lays its nine 1/d^2 bins between INCLUDE_RESOLUTION_RANGE's two values, not
between the extremes of the surviving data, and counts POSSIBLE against the
declared low limit - which is why its innermost shell reports the beam stop's
loss. Verified against a CORRECT.LP: all nine boundaries reproduce to the printed
precision from the declared 50 A, and not from the coarsest observed reflection.

Measured on stored merges of seven rotation datasets, small-molecule and protein,
re-scaled with --mode scale: the overall denominator moves by 0 to 2 reflections
out of 70,000-100,000, because on every one of them the coarsest reflection the
declared limit allows was itself measured - the corpus has no dataset whose stop
eats a whole low-resolution class. What does move is the shell grid: the
innermost boundary shifts by 0.1-0.4% in d (e.g. 7.21 -> 7.22 A), which changes
the innermost shell's counts by up to a few per cent and its R_meas by around
0.1 percentage points. Stored battery baselines for rmeas_lo must therefore be
regenerated, not compared across this commit.

Two decisions read merged completeness (the two-pass wrong-cell guard, which only
fires above 100.5% and only under -S); a larger denominator can only lower the
figure, so the guard can fire less often, never more.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
2026-09-02 16:39:40 +02:00

258 lines
12 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 "../image_analysis/scale_merge/HKLKey.h"
#include "../image_analysis/scale_merge/Merge.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.rlp = 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));
}
// --- 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);
}
}
Reference in New Issue View Git Blame Copy Permalink