From 3317eb6e88875dce0dccfdd256c63e04e640db52 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 7 Nov 2026 09:46:55 +0100 Subject: [PATCH] Correction surfaces: the resolution gauge counts observations The detector-frame surfaces are held to a mean of one within each resolution shell, because a factor that is a function of d alone is a null direction of their fit. That mean was weighted by 1/sigma^2, and sigma is smallest where little light arrives: on a detector part-shadowed by the pin the dark cells carry the shell's mean, and their large factors go onto the lit cells as small ones. How much of a ring is shadowed varies with resolution, so this bends the merged intensities along resolution: REFMAC's fitted scale against the model ran 1.48 / 0.80 / 1.16 in F (1.65 A / 2.7 A / inner) on a half-shadowed lysozyme sweep, with R as good as XDS's within every shell and 0.02 worse over all, and no shell-wise statistic can see it. Each observation now counts once. REFMAC R-work / R-free on the depositor's flags (d_min of the baseline), model as deposited, 31 sets chosen by the largest removed ramp and clamped cells, largest R vs per-shell R gap, plus controls: 5 better by more than 0.0015, 23 the same, 3 worse. half-shadowed lysozyme 0.3266 / 0.3306 -> 0.2996 / 0.3069 (XDS on the same model 0.3036 / 0.3083) soft X-ray insulin 0.2720 / 0.2841 -> 0.2423 / 0.2570 4p8r 0.2143 / 0.2518 -> 0.1873 / 0.2216 (depositor's data 0.2216) 6i3j -0.0029; insulin reference -0.0017; 2wnz +0.0025, 6jgh +0.0016, 5j23 +0.0015. The ramp removed on the first falls from 186x to 2x. The merge's own hash test records change as they must. A mean over the ring still lowers its lit part by the shadow it averages in (the new test reads 0.72 where 0.57 before and 1 would be exact); a median of the counts did not do better on the half-shadowed sweep (R-free 0.3090), the only one it was measured on. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQs733XDnbQQyUbV1vs --- .../scale_merge/RotationScaleMerge.cpp | 12 +++-- .../scale_merge/RotationScaleMergeGPU.cu | 2 +- .../scale_merge/RotationScaleMergeGPU.h | 6 +-- tests/RotationScaleMergeSegmentsTest.cpp | 53 +++++++++++++++++-- 4 files changed, 61 insertions(+), 12 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 79a02bd2c..8be21ff68 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4519,9 +4519,13 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int // that each shell still holds tens of thousands of observations. constexpr int NSHELL = 20; const int nshell = n_terms >= 100 * NSHELL ? NSHELL : 0; - // Weight of cell c in shell s: the statistical weight its observations carry there. Computed ONCE - // from the incoming sigmas rather than per round, so the projection is a fixed linear operator and - // the alternating fit below converges against it. + // Weight of cell c in shell s: the number of observations it holds there. Computed ONCE rather than + // per round, so the projection is a fixed linear operator and the alternating fit below converges + // against it. Not their statistical weight 1/sigma^2: sigma is smallest where little light arrives, + // so a cell behind a shadow would carry the shell's mean, its large factor going onto the lit cells + // as a small one. How much of a ring is shadowed varies with resolution, so that is a wave in the + // merged intensities, one no shell-wise statistic can see. Measured on a sweep with a third of the + // detector shadowed: 186x of ramp removed, 3x in between shells against the deposited model. std::vector shw(static_cast(ncell) * nshell, 0.0); std::vector shw_shell(nshell, 0.0), shw_cell(nshell > 0 ? ncell : 0, 0.0); std::vector group_shell(n_groups, 0); // the shell each ASU group sits in, for the gate @@ -4552,7 +4556,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int group_shell[t.group] = s; const double sc = static_cast(t.sigma) * t.corr; if (!(sc > 0.0)) continue; - shw[static_cast(t.cell) * nshell + s] += 1.0 / (sc * sc); + shw[static_cast(t.cell) * nshell + s] += 1.0; } } for (int c = 0; c < ncell && nshell > 0; ++c) diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 562b255d1..6288909d9 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -687,7 +687,7 @@ namespace { const auto &t = term[order[q]]; const double sc = __dmul_rn(double(t.sigma), double(t.corr)); if (!(sc > 0.0)) continue; - w = __dadd_rn(w, 1.0 / __dmul_rn(sc, sc)); + w = __dadd_rn(w, 1.0); } shw[b] = w; } diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index eb69648e2..210f17105 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -308,9 +308,9 @@ public: int block); // The resolution gauge's inputs for the terms SurfaceBuildTerms left: each group's shell (that of its - // last term in fulls order; length n_groups) and each (cell, shell) bin's statistical weight - // sum 1/(sigma corr)^2 in fulls order (length ncell * nshell), over nshell equal-occupancy shells in - // 1/d^2 - ApplyCellSurface's numbers to the bit. + // last term in fulls order; length n_groups) and each (cell, shell) bin's observation count + // (length ncell * nshell), over nshell equal-occupancy shells in 1/d^2 - ApplyCellSurface's numbers + // to the bit. void SurfaceShellWeights(int n_terms, int nshell, int32_t *group_shell, double *shw); // One subset of the terms (0 = even frames, 1 = odd, 2 = all), cut into the host's n_blocks diff --git a/tests/RotationScaleMergeSegmentsTest.cpp b/tests/RotationScaleMergeSegmentsTest.cpp index c5f20c7f9..1a1957398 100644 --- a/tests/RotationScaleMergeSegmentsTest.cpp +++ b/tests/RotationScaleMergeSegmentsTest.cpp @@ -202,8 +202,8 @@ TEST_CASE("RotationScaleMergeSegments_SegmentsMustCoverTheFrames", "[RotationSca TEST_CASE("RotationScaleMergeSegments_NoSegmentsIsTheEngineAsBefore", "[RotationScale]") { // The same synthetic sweep merged without segments gives, to the bit, what the engine gave before - // segments existed: FNV-1a over the merged I and sigma, recorded on the engine of rc181 e56b19abb - - // on the CPU, and on the GPU for the card model it was recorded on (the device path matches the host + // segments existed: FNV-1a over the merged I and sigma, recorded on the engine with the surface gauge + // counting observations - on the CPU, and on the GPU for the card model it was recorded on (the device path matches the host // path to the bit only where every reduction does; this one carries its own, so it is checked against // its own record). const std::vector sweeps{{40, 1.0, 0.0}, {40, 0.7, 3.0}}; @@ -226,9 +226,54 @@ TEST_CASE("RotationScaleMergeSegments_NoSegmentsIsTheEngineAsBefore", "[Rotation CHECK(r.segments.empty()); CHECK(r.merged.size() == 4106); if (!on_gpu) - CHECK(hash == 3376995616949734726ULL); + CHECK(hash == 15900277682318835699ULL); else if (get_gpu_description().find("RTX 5080") != std::string::npos) - CHECK(hash == 6594155396455850264ULL); + CHECK(hash == 6640579017825061097ULL); +} + +TEST_CASE("RotationScaleMerge_ShadowDoesNotBendTheResolutionShape", "[RotationScale]") { + // Every reflection lands on a ring of radius proportional to 1/d, at a random azimuth. Past a fixed + // detector column a shadow passes a fifth of the light and takes the background with it, so those + // observations have the small sigmas of a dark region - and it reaches only the rings wide enough to + // cross that column. The surface takes the shadow out; what it must not do is leave the lit part of + // the shadowed rings dimmer than the rings the shadow never touches. + const std::vector sweeps{{400, 1.0, 0.0}}; + auto data = MakeSweeps(sweeps, 21); + std::mt19937 rng(22); + std::uniform_real_distribution azimuth(0.0f, 2.0f * static_cast(M_PI)); + std::normal_distribution noise(0.0, 1.0); + for (auto &frame: data.frames) + for (auto &r: frame.reflections) { + const float radius = 300.0f / r.d; + const float phi = azimuth(rng); + r.predicted_x = 500.0f + radius * std::cos(phi); + r.predicted_y = 500.0f + radius * std::sin(phi); + if (r.predicted_x > 590.0f) { + const double expected = data.I_true.at(FriedelKey(r.h, r.k, r.l)) * 0.2; + const double sigma = std::sqrt(expected + 1.0); + r.I = static_cast(expected + sigma * noise(rng)); + r.sigma = static_cast(sigma); + r.var_bkg = static_cast(sigma * sigma - std::max(0.0, static_cast(r.I))); + r.bkg = 1.0f; + } + } + const auto result = Merge(data, SyntheticExperiment(true), {}); + + // / in the rings the shadow misses (d > 3.4 A) against those it crosses + double in_merged = 0.0, in_true = 0.0, out_merged = 0.0, out_true = 0.0; + for (const auto &m: result.merged) { + const auto t = data.I_true.find(FriedelKey(m.h, m.k, m.l)); + REQUIRE(t != data.I_true.end()); + (m.d > 3.4 ? in_merged : out_merged) += m.I; + (m.d > 3.4 ? in_true : out_true) += t->second; + } + REQUIRE(in_true > 0.0); + REQUIRE(out_true > 0.0); + const double shape = (out_merged / out_true) / (in_merged / in_true); + // A mean over the ring still lowers its lit part by the shadow it averages in: 0.72 with every + // observation counted once, 0.57 where the dark ones' small sigmas gave them the weight. + CHECK(shape > 0.65); + CHECK(shape < 1.25); } // A sweep in which no rocking event is captured whole enough to make a full - as in one exposed 0.2 deg out