RotationScaleMerge: fulls-only per-frame scaling (prototype), pooled over sparse frames

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

common/ScalingSettings.cpp

@@ -208,6 +208,15 @@ double ScalingSettings::GetSmoothGDegrees() const {
     return smooth_g_deg;
 }
 
+ScalingSettings &ScalingSettings::ScalePartials(bool input) {
+    scale_partials = input;
+    return *this;
+}
+
+bool ScalingSettings::GetScalePartials() const {
+    return scale_partials;
+}
+
 ScalingSettings &ScalingSettings::RelativeBDegrees(double input) {
     if (input < 0)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Relative-B batch width must be non-negative");

common/ScalingSettings.h

@@ -98,6 +98,7 @@ class ScalingSettings {
     // degrees (like XDS DELPHI), converted to an odd frame window from the oscillation step; this keeps
     // the smoothing physical (independent of frame slicing). 0 = off. A no-op without rot3d.
     double smooth_g_deg = 0.0;
+    bool scale_partials = true;
 
     // Per-batch relative-B on the rot3d fulls (beyond the single global decay slope): bin frames into
     // rotation-range batches of this width in degrees and refine one relative Debye-Waller B per batch, so
@@ -143,6 +144,7 @@ public:
     ScalingSettings& IceMinScore(float input);
     ScalingSettings& IceMinSpotRatio(float input);
     ScalingSettings& SmoothGDegrees(double input);
+    ScalingSettings& ScalePartials(bool input);
     ScalingSettings& RelativeBDegrees(double input);
 
     ScalingSettings& RfreeFraction(double input);
@@ -182,6 +184,7 @@ public:
     [[nodiscard]] float GetIceMinScore() const;
     [[nodiscard]] float GetIceMinSpotRatio() const;
     [[nodiscard]] double GetSmoothGDegrees() const;
+    [[nodiscard]] bool GetScalePartials() const;
     [[nodiscard]] double GetRelativeBDegrees() const;
 
     [[nodiscard]] double GetRfreeFraction() const;

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -1407,17 +1407,34 @@ void RotationScaleMerge::ReduceScalingGroupMeans(int n_groups, const std::vector
     });
 }
 
+int RotationScaleMerge::PoolHalfWidth(const std::vector<int32_t> &fcount, int f) {
+    const int n = static_cast<int>(fcount.size());
+    int64_t total = fcount[f];
+    int h = 0;
+    while (total < static_cast<int64_t>(MIN_REFLECTIONS) && h < n) {
+        ++h;
+        if (f - h >= 0) total += fcount[f - h];
+        if (f + h < n) total += fcount[f + h];
+    }
+    return h;
+}
+
 template <class O>
 void RotationScaleMerge::FitPerFrameG(const std::vector<O> &obs, const std::vector<int32_t> &fstart,
                                       const std::vector<int32_t> &fcount,
                                       const std::vector<double> &group_mean_in,
                                       bool unity, std::vector<double> &g) {
     std::vector<uint8_t> scaled(fstart.size(), 0);
-    ParallelFor(static_cast<int>(fstart.size()), nthreads, [&](int f) {
+    const int n_fr = static_cast<int>(fstart.size());
+    ParallelFor(n_fr, nthreads, [&](int f) {
+        // A full is a whole rocking event filed under one frame, so on a sparse sweep a frame holds a
+        // handful of them - too few to fit a scale on, and the frame was left unscaled, which on such a
+        // sweep was every frame. Those frames are fitted over the nearest frames on either side that
+        // together hold enough; a frame with enough of its own is fitted on its own, as before.
+        const int h = unity && pool_sparse_fulls ? PoolHalfWidth(fcount, f) : 0;
         std::vector<ScaleObs> so;
-        so.reserve(fcount[f]);
-        const int lo = fstart[f], hi = fstart[f] + fcount[f];
-        for (int i = lo; i < hi; ++i) {
+        for (int j = std::max(0, f - h); j <= std::min(n_fr - 1, f + h); ++j)
+        for (int i = fstart[j]; i < fstart[j] + fcount[j]; ++i) {
             const auto &o = obs[i];
             if (o.group < 0) continue;
             if (o.on_ice) continue;
@@ -5607,9 +5624,15 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
                     "curves span {} frames)", smooth_window, smooth_window * osc_deg, smooth_g_deg,
                     rocking_event_frames_at_start);
     }
+    // Without partial scaling every frame keeps G = 1 here and its scale comes from the fulls alone
+    // (scale-fulls below), as in XDS: a partial's scale and an error of the partiality model are the
+    // same thing within a rocking curve, and on a sparse fine-sliced sweep the fit takes the one for
+    // the other.
+    const bool scale_partials = s.GetScalePartials();
+    pool_sparse_fulls = !scale_partials;
     ScalingLoopOutcome partial_loop;
 #ifdef JFJOCH_USE_CUDA
-    if (gpu_active_) {
+    if (gpu_active_ && scale_partials) {
         // The scaling loop runs on the GPU one iteration per call, and corr stays RESIDENT across
         // scaling -> smooth-G -> CC -> combine (and across passes, exactly as the old host round-trip
         // did). Only the per-frame G/scaled come back after each iteration, for the gauge pin and the
@@ -5630,7 +5653,11 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
         scaled_on_gpu = true;
     }
 #endif
-    if (!scaled_on_gpu) {
+    if (!scale_partials) {
+        std::fill(g_partial.begin(), g_partial.end(), 1.0);
+        frame_scaled_scratch.assign(n_frames, 1);
+        partial_loop.converged = true;
+    } else if (!scaled_on_gpu) {
         const PartialLoopKey key{x.GetSpaceGroupOrP1().xhm(), merge_friedel, d_min_limit, d_max_limit,
                                  min_partiality, smooth_window, scaling_iter};
         const auto memo = std::find_if(partial_loop_memos.begin(), partial_loop_memos.end(),
@@ -5920,7 +5947,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
             std::vector<uint8_t> scaled_dev(n_frames);
             fulls_loop = RunScalingLoop("fulls", g_full, f_count, smooth_window,
                                           /*release_restraint=*/true, [&] {
-                gpu_->ScaleFulls(1, min_partiality);
+                gpu_->ScaleFulls(1, min_partiality, pool_sparse_fulls);
                 gpu_->GetG(g_dev.data(), scaled_dev.data());
                 for (int f = 0; f < n_frames; ++f)
                     if (scaled_dev[f]) g_full[f] = g_dev[f];

image_analysis/scale_merge/RotationScaleMerge.h

@@ -43,6 +43,10 @@
 // Stills use the per-image ScaleOnTheFly (fixed partiality) instead.
 class RotationScaleMerge {
 public:
+    // How many frames on either side a frame's fulls scale is fitted over (used without partial
+    // scaling): none when the frame holds MIN_REFLECTIONS fulls itself, else the fewest that bring the
+    // pooled count there.
+    static int PoolHalfWidth(const std::vector<int32_t> &fcount, int f);
     struct Result {
         std::vector<MergedReflection> merged;
         MergeStatistics statistics;
@@ -213,6 +217,9 @@ private:
     // this is how far they may go before giving up and saying so.
     int scaling_iter = 100;
     bool scale_fulls = true;
+    // Without partial scaling the fulls carry the whole per-frame scale, and a sparse frame is fitted
+    // over its neighbours (PoolHalfWidth).
+    bool pool_sparse_fulls = false;
     bool refine_decay_b = false;   // per-time-block Debye-Waller decay correction (radiation damage)
     int  absorption_iter = 0;      // >0: fit a goniometer-frame absorption surface over this many iterations
     int  modulation_iter = 0;      // >0: fit a detector-plane modulation (flat-field) surface, this many iterations

image_analysis/scale_merge/RotationScaleMergeGPU.cu

@@ -123,16 +123,31 @@ namespace {
                                        const int32_t *__restrict__ frame_count,
                                        const float *__restrict__ I, const double *__restrict__ inv_sigma,
                                        const float *__restrict__ sco_coeff, const uint8_t *__restrict__ sco_ok,
-                                       const int32_t *__restrict__ perm,
+                                       const int32_t *__restrict__ perm, bool pool,
                                        double *__restrict__ g, uint8_t *__restrict__ scaled) {
         const int f = blockIdx.x;
         if (f >= n_frames) return;
-        const int lo = frame_start[f], hi = frame_start[f] + frame_count[f];
+        // pool (the fulls): a frame holding fewer than MIN_REFLECTIONS is fitted over the fewest frames
+        // on either side that bring the count there - RotationScaleMerge::PoolHalfWidth, same rule.
+        __shared__ int s_h;
+        if (threadIdx.x == 0) {
+            long total = frame_count[f];
+            int h = 0;
+            while (pool && total < long(MIN_REFLECTIONS) && h < n_frames) {
+                ++h;
+                if (f - h >= 0) total += frame_count[f - h];
+                if (f + h < n_frames) total += frame_count[f + h];
+            }
+            s_h = h;
+        }
+        __syncthreads();
+        const int j_lo = max(0, f - s_h), j_hi = min(n_frames - 1, f + s_h);
         __shared__ double sh[BLK];
 
         long cnt_local = 0;
-        for (int i = lo + threadIdx.x; i < hi; i += blockDim.x)
-            if (sco_ok[perm ? perm[i] : i]) ++cnt_local;
+        for (int j = j_lo; j <= j_hi; ++j)
+            for (int i = frame_start[j] + threadIdx.x; i < frame_start[j] + frame_count[j]; i += blockDim.x)
+                if (sco_ok[perm ? perm[i] : i]) ++cnt_local;
         const double cnt = BlockReduceSum(double(cnt_local), sh);
         __shared__ double s_cnt;
         if (threadIdx.x == 0) s_cnt = cnt;
@@ -140,15 +155,16 @@ namespace {
         if (s_cnt < MIN_REFLECTIONS) return;   // leave g[f]/scaled[f] as-is
 
         double num = 0.0, den = 0.0;
-        for (int i = lo + threadIdx.x; i < hi; i += blockDim.x) {
-            const int a = perm ? perm[i] : i;
-            if (!sco_ok[a]) continue;
-            const double coeff = sco_coeff[a];
-            const double w = inv_sigma[a];
-            const double w2 = w * w;
-            num += w2 * coeff * double(I[a]);
-            den += w2 * coeff * coeff;
-        }
+        for (int j = j_lo; j <= j_hi; ++j)
+            for (int i = frame_start[j] + threadIdx.x; i < frame_start[j] + frame_count[j]; i += blockDim.x) {
+                const int a = perm ? perm[i] : i;
+                if (!sco_ok[a]) continue;
+                const double coeff = sco_coeff[a];
+                const double w = inv_sigma[a];
+                const double w2 = w * w;
+                num += w2 * coeff * double(I[a]);
+                den += w2 * coeff * coeff;
+            }
         const double tnum = BlockReduceSum(num, sh); __syncthreads();
         const double tden = BlockReduceSum(den, sh);
         if (threadIdx.x == 0) {
@@ -1081,7 +1097,7 @@ void RotationScaleMergeGPU::ScalePartials(int iters, double min_partiality, bool
             d.sco_coeff.get(), d.sco_ok.get());
         CudaCheck(cudaGetLastError(), "PrepScaleObsKernel launch");
         FitPerFrameGKernel<<<d.n_frames, BLK, 0, impl_->s()>>>(d.n_frames, d.frame_start.get(), d.frame_count.get(),
-            d.I.get(), d.inv_sigma.get(), d.sco_coeff.get(), d.sco_ok.get(), nullptr, d.g.get(), d.scaled.get());
+            d.I.get(), d.inv_sigma.get(), d.sco_coeff.get(), d.sco_ok.get(), nullptr, /*pool=*/false, d.g.get(), d.scaled.get());
         CudaCheck(cudaGetLastError(), "FitPerFrameGKernel launch");
         UpdateCorrKernel<<<upd_blocks, BLK, 0, impl_->s()>>>(d.n_obs, d.frame.get(), d.prescaling_corr.get(), d.partiality.get(),
             d.g.get(), d.scaled.get(), d.corr.get());
@@ -1497,7 +1513,7 @@ void RotationScaleMergeGPU::ResetFullsScale() {
     CudaCheck(cudaStreamSynchronize(impl_->s()), "reset fulls scale sync");
 }
 
-void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality) {
+void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality, bool pool) {
     DeviceGuard guard(impl_->device, impl_->available);
     auto &d = *impl_;
     const int nf = d.n_fulls;
@@ -1523,7 +1539,7 @@ void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality) {
         CudaCheck(cudaGetLastError(), "PrepScaleObsKernel launch");
         FitPerFrameGKernel<<<d.n_frames, BLK, 0, impl_->s()>>>(d.n_frames,
             d.f_frame_start.get(), d.f_frame_count.get(), d.f_I.get(), d.f_inv_sigma.get(),
-            d.f_sco_coeff.get(), d.f_sco_ok.get(), d.f_frame_perm.get(), d.g.get(), d.scaled.get());
+            d.f_sco_coeff.get(), d.f_sco_ok.get(), d.f_frame_perm.get(), pool, d.g.get(), d.scaled.get());
         CudaCheck(cudaGetLastError(), "FitPerFrameGKernel launch");
         UpdateCorrKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(nf, d.f_frame.get(), d.f_rlp.get(), d.f_partiality.get(),
             d.g.get(), d.scaled.get(), d.f_corr.get());

image_analysis/scale_merge/RotationScaleMergeGPU.h

@@ -184,7 +184,9 @@ public:
     // Run `iters` of the Unity scaling loop on the resident fulls (reduce group means -> per-frame LS G
     // -> update corr), in place on the fulls' working corr. Requires SetFullsFrameCSR + SetFullsGroups
     // and ResetFullsScale.
-    void ScaleFulls(int iters, double min_partiality);
+    // pool: a frame with fewer than MIN_REFLECTIONS fulls is fitted over its neighbours
+    // (RotationScaleMerge::PoolHalfWidth).
+    void ScaleFulls(int iters, double min_partiality, bool pool);
 
     // The fulls' counterpart of SmoothCorr: f_corr[i] *= ratio[f_frame[i]] where apply[f].
     void SmoothFullsCorr(const uint8_t *apply, const double *ratio);

rugnux/rugnux_cli.cpp

@@ -157,6 +157,7 @@ void print_usage() {
     std::cout << "   --write-process-h5              Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl;
     std::cout << "   --finalist-ledger               Report the full-resolution evidence for each space group the search considered, not only the one it adopted (report-only; the decision is unchanged)" << std::endl;
     std::cout << "   --developer                     Write the full <prefix>_report.txt: the pipeline-internal keys (the anisotropy gate, the space-group candidate and operator tables, the model-fit null, the sweep internals) and the long explanations, which the default report leaves out. Nothing is computed differently - the same report, rendered in full" << std::endl;
+    std::cout << "   --no-scale-partials             rot3d: do not fit per-frame scales on the partials; the per-frame scale comes from the fulls (scale-fulls) alone, as in XDS" << std::endl;
     std::cout << "   --smooth-g[=deg]                rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl;
     std::cout << "   --relative-b[=deg]              rot3d: fit a per-batch relative-B (beyond the single decay slope) over deg-degree batches; cross-validated (default: 10 deg when bare; off otherwise)" << std::endl;
     std::cout << "   --no-scaling-corrections        rot3d: disable the (default-on) decay + absorption + modulation correction surfaces fitted on the fulls after scale-fulls" << std::endl;
@@ -282,6 +283,7 @@ enum {
     OPT_MOSAICITY,
     OPT_PREDICTION_MOSAICITY,
     OPT_SMOOTH_G,
+    OPT_NO_SCALE_PARTIALS,
     OPT_RELATIVE_B,
     OPT_NO_SCALING_CORRECTIONS,
     OPT_NO_EXPECTED_VARIANCE_MERGE,
@@ -351,6 +353,7 @@ static option long_options[] = {
     {"finalist-ledger", no_argument, nullptr, OPT_FINALIST_LEDGER},
     {"developer", no_argument, nullptr, OPT_DEVELOPER},
     {"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G},
+    {"no-scale-partials", no_argument, nullptr, OPT_NO_SCALE_PARTIALS},
     {"relative-b", optional_argument, nullptr, OPT_RELATIVE_B},
     {"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS},
     {"no-expected-variance-merge", no_argument, nullptr, OPT_NO_EXPECTED_VARIANCE_MERGE},
@@ -773,6 +776,7 @@ static int RunRugnux(int argc, char **argv) {
     std::optional<float> beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, rot3_rad, polarization_factor;
     bool detector_mirror_y = false;
     int64_t detector_quarter_turns = 0;
+    bool no_scale_partials = false;           // --no-scale-partials (prototype)
     std::optional<double> smooth_g_deg_arg;   // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise
     std::optional<double> relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given
     bool no_scaling_corrections = false;   // --no-scaling-corrections: disable rot3d decay+absorption+modulation surfaces
@@ -1207,6 +1211,9 @@ static int RunRugnux(int argc, char **argv) {
         case OPT_DEVELOPER:
             provenance.developer = true;
             break;
+        case OPT_NO_SCALE_PARTIALS:
+            no_scale_partials = true;
+            break;
         case OPT_SMOOTH_G:
             smooth_g_deg_arg = optarg ? parse_double_arg(optarg, "--smooth-g", logger) : SMOOTH_G_DEFAULT_DEG;
             break;
@@ -1652,6 +1659,7 @@ static int RunRugnux(int argc, char **argv) {
             outlier_reject_nsigma.value_or(scaling_settings.GetOutlierRejectNsigma()));
         scaling_settings.ScaleFulls(scale_fulls_arg.value_or(scaling_settings.GetScaleFulls()));
         scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees()));
+    if (no_scale_partials) scaling_settings.ScalePartials(false);
         scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0));   // opt-in only; default off
         if (no_scaling_corrections)
             scaling_settings.CorrectionSurfaces(false);
@@ -2538,6 +2546,7 @@ static int RunRugnux(int argc, char **argv) {
     ScalingSettings scaling_settings = RugnuxDefaultScalingSettings(rotation_indexing);
     scaling_settings.ScaleFulls(scale_fulls);
     scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees()));
+    if (no_scale_partials) scaling_settings.ScalePartials(false);
     scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0));   // opt-in only; default off
     if (no_scaling_corrections)
         scaling_settings.CorrectionSurfaces(false);

tests/RotationScaleWalkTest.cpp

@@ -5,6 +5,7 @@
 #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,
@@ -92,3 +93,15 @@ TEST_CASE("WalkRotationScale_StoredAnglesStand", "[RotationScale]") {
         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);
+}