From f69339ce6387234c8e3fdde14f9918238dd50725 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 23:08:12 +0200 Subject: [PATCH] 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) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- common/ScalingSettings.cpp | 9 ++++ common/ScalingSettings.h | 3 ++ .../scale_merge/RotationScaleMerge.cpp | 41 +++++++++++++--- .../scale_merge/RotationScaleMerge.h | 7 +++ .../scale_merge/RotationScaleMergeGPU.cu | 48 ++++++++++++------- .../scale_merge/RotationScaleMergeGPU.h | 4 +- rugnux/rugnux_cli.cpp | 9 ++++ tests/RotationScaleWalkTest.cpp | 13 +++++ 8 files changed, 110 insertions(+), 24 deletions(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 2016afba5..4fabbd44d 100644 --- a/common/ScalingSettings.cpp +++ b/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"); diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 198c6366c..e9b2cd72e 100644 --- a/common/ScalingSettings.h +++ b/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; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 81dabcd22..2b71049b4 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/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 &fcount, int f) { + const int n = static_cast(fcount.size()); + int64_t total = fcount[f]; + int h = 0; + while (total < static_cast(MIN_REFLECTIONS) && h < n) { + ++h; + if (f - h >= 0) total += fcount[f - h]; + if (f + h < n) total += fcount[f + h]; + } + return h; +} + template void RotationScaleMerge::FitPerFrameG(const std::vector &obs, const std::vector &fstart, const std::vector &fcount, const std::vector &group_mean_in, bool unity, std::vector &g) { std::vector scaled(fstart.size(), 0); - ParallelFor(static_cast(fstart.size()), nthreads, [&](int f) { + const int n_fr = static_cast(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 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 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]; diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 79278560f..30da42438 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/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 &fcount, int f); struct Result { std::vector 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 diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 66aeabacb..06b85d5bd 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/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<<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<<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<<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<<s()>>>(nf, d.f_frame.get(), d.f_rlp.get(), d.f_partiality.get(), d.g.get(), d.scaled.get(), d.f_corr.get()); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 766e2f640..98aa1b18a 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/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); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 4eb457f6b..9cb7403f5 100644 --- a/rugnux/rugnux_cli.cpp +++ b/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 _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 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 smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise std::optional 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); diff --git a/tests/RotationScaleWalkTest.cpp b/tests/RotationScaleWalkTest.cpp index 562d46a24..f01334b3f 100644 --- a/tests/RotationScaleWalkTest.cpp +++ b/tests/RotationScaleWalkTest.cpp @@ -5,6 +5,7 @@ #include #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 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); +}