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); +}