Prototype switch --no-scale-partials: skip the per-frame scaling of the partials and let the fulls (scale-fulls) carry the per-frame scale, as XDS does. Without partial scaling a frame whose fulls number fewer than MIN_REFLECTIONS is fitted over the nearest frames on either side that together hold enough (PoolHalfWidth), on the host and in the GPU kernel. Default path unchanged (p.mtz md5 identical on the three profiling sets). Why: on fine-sliced small-molecule sweeps the per-frame partial scale and the partiality model are degenerate within a rocking curve, and the fit swings G 0.23..1.2 with a 180 deg period (XDS's own frame scale: 0.79..0.99). That imprints an hkl-dependent bias common to all equivalents, which R_meas/CC1/2/ISa cannot see but a refinement against the known structure does. And scale-fulls never fitted a frame on such data: a full is filed under one frame, about 8 per frame, below MIN_REFLECTIONS, so every frame kept G = 1. Measured with SHELXL refining the COD structures (R1 >4sig), default -> switch: aspirin 20 keV 0.096 -> 0.062, aspirin 25 keV 0.094 -> 0.046, citric acid 0.161 -> 0.109, HEPES 0.090 -> 0.070 (XDS 0.030-0.038). Proteins lose ISa with the switch (myob 9.1 -> 7.6, cytc 25.8 -> 13.7), so it is not a default; the choice is to be made from the data. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB
108 lines
4.8 KiB
C++
108 lines
4.8 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <cmath>
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include "../rugnux/Rugnux.h"
|
|
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
|
|
|
|
namespace {
|
|
// A synthetic sweep whose stage turned `true_scale` times the stored angles. Scored at a scale k,
|
|
// the validation spots stay on the lattice as long as the angles track the rotation, and the share
|
|
// that does falls off with the relative rate error; a wrong spindle angle keeps half a percent.
|
|
struct SyntheticSweep {
|
|
double true_scale;
|
|
int64_t spots = 35000;
|
|
int index_calls = 0;
|
|
int refit_calls = 0;
|
|
|
|
ValidationSpotEvidence IndexAt(float k) {
|
|
++index_calls;
|
|
const double error = std::fabs(true_scale / k - 1.0);
|
|
const double on = 0.9 * std::max(0.0, 1.0 - 30.0 * error);
|
|
return ValidationSpotEvidence{spots, std::llround(on * spots), std::llround(0.005 * spots)};
|
|
}
|
|
|
|
// The post-refinement at k, relative to k. It reads only 70 % of the error that is left: it
|
|
// sees only the frames the angles at k still track.
|
|
std::optional<double> RefitAt(float k) {
|
|
++refit_calls;
|
|
return 1.0 + 0.7 * (true_scale / k - 1.0);
|
|
}
|
|
|
|
RotationScaleWalk Walk(double first_fit) {
|
|
return WalkRotationScale(first_fit, [this](float k) { return IndexAt(k); },
|
|
[this](float k) { return RefitAt(k); }, 8);
|
|
}
|
|
};
|
|
}
|
|
|
|
TEST_CASE("ValidationEvidencePrefers", "[RotationScale]") {
|
|
const ValidationSpotEvidence base{10000, 3000, 50};
|
|
// 1 % more of the spots beyond chance is under the noise of two 30 % shares over 10000 spots
|
|
// (sqrt(2 * 0.3 * 0.7 / 10000) = 0.65 %, times 3.29); 5 % is well over it.
|
|
CHECK_FALSE(ValidationEvidencePrefers(base, {10000, 3100, 50}));
|
|
CHECK(ValidationEvidencePrefers(base, {10000, 3500, 50}));
|
|
// A candidate is judged against its own null: more spots on the lattice bought by a null that
|
|
// rose just as much is no gain.
|
|
CHECK_FALSE(ValidationEvidencePrefers(base, {10000, 3500, 550}));
|
|
// Never against itself, and nothing that scored nothing wins.
|
|
CHECK_FALSE(ValidationEvidencePrefers(base, base));
|
|
CHECK_FALSE(ValidationEvidencePrefers(base, {}));
|
|
CHECK(ValidationEvidencePrefers({}, base));
|
|
}
|
|
|
|
TEST_CASE("WalkRotationScale_ReachesTheFixedPoint", "[RotationScale]") {
|
|
// A stage 3 % slow. The first fit reads 70 % of that; each refit at the adopted scale reads 70 %
|
|
// of what is left, and the walk goes on as long as the validation frames prefer the new scale.
|
|
SyntheticSweep sweep{0.97};
|
|
const auto walk = sweep.Walk(sweep.RefitAt(1.0f).value());
|
|
CHECK(walk.scale == Catch::Approx(0.97).margin(0.001));
|
|
CHECK(walk.scale != 1.0f);
|
|
CHECK(walk.evidence.on_lattice > sweep.IndexAt(1.0f).on_lattice);
|
|
CHECK(sweep.refit_calls > 2);
|
|
CHECK_FALSE(walk.trail.empty());
|
|
}
|
|
|
|
TEST_CASE("WalkRotationScale_StoredAnglesStand", "[RotationScale]") {
|
|
SECTION("A healthy stage: a fit off by noise scores no better than the stored angles") {
|
|
SyntheticSweep sweep{1.0};
|
|
const auto walk = sweep.Walk(1.0002);
|
|
CHECK(walk.scale == 1.0f);
|
|
CHECK(sweep.index_calls == 2); // the stored angles and the fit, nothing more
|
|
CHECK(sweep.refit_calls == 0);
|
|
}
|
|
SECTION("A real but small error the spots cannot resolve beyond their noise") {
|
|
SyntheticSweep sweep{0.999};
|
|
sweep.spots = 400;
|
|
const auto walk = sweep.Walk(0.9993);
|
|
CHECK(walk.scale == 1.0f);
|
|
}
|
|
SECTION("A fit that tracks something other than the rotation scores worse, and is refused") {
|
|
SyntheticSweep sweep{1.0};
|
|
const auto walk = sweep.Walk(0.98);
|
|
CHECK(walk.scale == 1.0f);
|
|
CHECK(sweep.refit_calls == 0);
|
|
}
|
|
SECTION("A fit of exactly one asks for no probe at all") {
|
|
SyntheticSweep sweep{1.0};
|
|
const auto walk = sweep.Walk(1.0);
|
|
CHECK(walk.scale == 1.0f);
|
|
CHECK(walk.trail.empty());
|
|
CHECK(sweep.index_calls == 0);
|
|
}
|
|
}
|
|
|
|
TEST_CASE("PoolHalfWidth_FitsSparseFramesOverTheirNeighbours", "[RotationScale]") {
|
|
// 20 fulls of its own: fitted on its own.
|
|
CHECK(RotationScaleMerge::PoolHalfWidth({0, 20, 0}, 1) == 0);
|
|
// 5 per frame: two frames either side bring 25 >= 20.
|
|
const std::vector<int32_t> five(11, 5);
|
|
CHECK(RotationScaleMerge::PoolHalfWidth(five, 5) == 2);
|
|
// At the end of the sweep the window grows on the one side there is.
|
|
CHECK(RotationScaleMerge::PoolHalfWidth(five, 0) == 3);
|
|
// A sweep that never holds enough stops at its length.
|
|
CHECK(RotationScaleMerge::PoolHalfWidth({1, 1, 1}, 1) == 3);
|
|
}
|