Rotation scaling: penalised per-frame scale of the fulls; fulls-only partials get flux and partiality
Three defects in the per-frame scaling of sparse (small-molecule, weak) rotation sweeps: 1. The fulls' per-frame scale pooled sparse frames by their RAW full count, but in a high-symmetry group with many systematic absences most usable counts stayed under the minimum, so most frames were never fitted and kept corr = 1 beside pinned, fitted frames (two gauges in one reference). The release step then read the gauge offset (a constant 134x on a cubic Ia-3d small-molecule sweep) as signal and gave each frame exp(kept_f * 4.9) - the e^-5 errors on a quarter of the frames, R1 0.62. A fixed box window also cannot follow a 100x absorption ramp over a few degrees, and frames under the credible floor, exempt from the window, ran away to 1e-7. Now each round fits every frame on its own fulls (no pooling, no minimum), and the scale is a penalised second-difference smoother of log G (Whittaker/Eilers), each frame at the information of its fit, lambda by cross-validation over blocks one rocking curve wide (interleaved single frames leak through shared rocking curves and chose to follow every frame). After convergence the existing ShrinkToRestrained hands back the per-frame deviation its neighbour shares. Pooling and the box window are gone from the fulls loop; the partials loop is unchanged. 2. With partial scaling off (< 50 rocking events per frame) the partials kept the integration-time corr: no incident-flux correction and not the partiality of the ingest-smoothed geometry, because only the partial scaling loop rewrote corr. They now get corr = prescaling_corr / partiality at G = 1 (host, and a device kernel). 3. The flux meter (per-frame mean background) jumped 30x between neighbouring frames of a sparse sweep - on a few reflections it measures which reflections the frame holds. It is read through the same smoother at the precision of each frame's mean. SHELXL R1(>4sig) against the published structures, rc174-cand -> this, on the in-house small-molecule sweeps: cubic Ia-3d 0.615 -> 0.119 (XDS 0.088); four organic sweeps (monoclinic / orthorhombic) 0.0515 -> 0.0459, 0.0451 -> 0.0396, 0.0747 -> 0.0709, 0.0452 -> 0.0413; ISa up to 11.6 -> 31. Raw flux instead of smoothed costs 0.002-0.003 R1 on the first two. A weak, decaying protein sweep on the fulls-only path: ISa 6.7 -> 27.1. CPU and GPU paths agree. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB
This commit is contained in:
@@ -94,14 +94,43 @@ TEST_CASE("WalkRotationScale_StoredAnglesStand", "[RotationScale]") {
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
TEST_CASE("SmoothLogScale_FollowsInformationBridgesGaps", "[RotationScale]") {
|
||||
const int n = 60;
|
||||
// A ramp is no curvature, so any amount of smoothing keeps it exactly.
|
||||
std::vector<double> ramp(n), J(n, 1.0);
|
||||
for (int f = 0; f < n; ++f) ramp[f] = -0.1 * f;
|
||||
auto x = RotationScaleMerge::SmoothLogScale(ramp, J, 1e6);
|
||||
for (int f = 0; f < n; ++f) CHECK(x[f] == Catch::Approx(ramp[f]).margin(1e-6));
|
||||
// A stretch with no information is bridged by the straight line through its neighbours.
|
||||
std::vector<double> Jgap(J);
|
||||
for (int f = 20; f < 40; ++f) Jgap[f] = 0.0;
|
||||
x = RotationScaleMerge::SmoothLogScale(ramp, Jgap, 1.0);
|
||||
CHECK(x[30] == Catch::Approx(-3.0).margin(1e-6));
|
||||
// A frame with far more information than its neighbours keeps its own value.
|
||||
std::vector<double> y(n, 0.0), Jone(n, 1.0);
|
||||
y[30] = 1.0; Jone[30] = 1e6;
|
||||
x = RotationScaleMerge::SmoothLogScale(y, Jone, 10.0);
|
||||
CHECK(x[30] == Catch::Approx(1.0).margin(1e-3));
|
||||
// Fewer than two frames with information: nothing to smooth against.
|
||||
std::vector<double> Jsingle(n, 0.0);
|
||||
Jsingle[5] = 1.0;
|
||||
CHECK(RotationScaleMerge::SmoothLogScale(y, Jsingle, 1.0) == y);
|
||||
}
|
||||
|
||||
TEST_CASE("ChooseLogScaleSmoothing_SmoothsNoiseFollowsSignal", "[RotationScale]") {
|
||||
const int n = 400;
|
||||
std::vector<double> J(n, 1.0), noisy(n), step(n);
|
||||
// Deterministic noise about a flat scale: the chosen curve is close to flat.
|
||||
for (int f = 0; f < n; ++f) noisy[f] = 0.3 * std::sin(12.9898 * f) * std::cos(78.233 * f);
|
||||
const double l_noise = RotationScaleMerge::ChooseLogScaleSmoothing(noisy, J, 1);
|
||||
const auto flat = RotationScaleMerge::SmoothLogScale(noisy, J, l_noise);
|
||||
double rms = 0.0;
|
||||
for (double v : flat) rms += v * v;
|
||||
CHECK(std::sqrt(rms / n) < 0.05);
|
||||
// A precise slow wave is followed.
|
||||
for (int f = 0; f < n; ++f) step[f] = 2.0 * std::sin(f / 30.0);
|
||||
const double l_wave = RotationScaleMerge::ChooseLogScaleSmoothing(step, J, 1);
|
||||
const auto wave = RotationScaleMerge::SmoothLogScale(step, J, l_wave);
|
||||
CHECK(wave[47] == Catch::Approx(step[47]).margin(0.01));
|
||||
CHECK(l_wave < l_noise);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user