From 2c928b27cd615f85a3c22cedd18a46c467139f44 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 3 Jul 2026 07:27:33 +0200 Subject: [PATCH] RotationScaleMerge: GPU 3D-combine (CUDA port, phase 2 step 1) Port Combine() (partials->fulls) to CUDA, mirroring process_rawrun bit-for-bit: one thread per raw-hkl run splits its usable partials into rocking events (frame gap <= 2), pools background, seeds F, runs 3 de-biased Poisson reweights and the capture-uncertainty term. Emission is deterministic - a count pass, a host exclusive prefix sum for per-run offsets, then an emit pass at those offsets - so fulls come out in raw-run-major/event order, identical to the CPU path; both pass instantiations share the same arithmetic so count == emit exactly. Dmax/Dmin/Fmax reproduce std::max/min NaN semantics (not fmax) for parity. Validated across the 18-crystal rotation battery: all 15 deterministic crystals (P1/P2/C2/H3/I23/P41212/P222/P422) match the CPU combine exactly on SG/ISa/CC1.2/ completeness and run-to-run (fulls count bit-identical); the 3 upstream-nondet crystals vary from GPU-prediction overflow, not the combine. Gated opt-in behind JFJOCH_RSM_GPU_COMBINE (default = CPU combine): combine alone is timing-neutral because the shared 1.2M SortFullsByFrame std::sort dominates and the fulls round-trip adds a copy - it only pays off once the fulls stay resident for scale-fulls + merge. Also add JFJOCH_RSM_NO_GPU master switch to force the CPU fallback (incl. phase-1 scaling) from one binary for A/B parity. SortFullsByFrame extracted from the Combine tail and shared by both paths. Co-Authored-By: Claude Opus 4.8 --- .../scale_merge/RotationScaleMerge.cpp | 52 +++- .../scale_merge/RotationScaleMerge.h | 11 +- .../scale_merge/RotationScaleMergeGPU.cu | 270 ++++++++++++++++++ .../scale_merge/RotationScaleMergeGPU.h | 26 ++ 4 files changed, 353 insertions(+), 6 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 5eec365e..4dd54dad 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -243,20 +244,27 @@ void RotationScaleMerge::Ingest() { // Bring the partial-scaling loop onto the GPU when one is present. Upload the immutable per-obs // fields once (corr lives on the device, refreshed each pass); the CPU keeps the sort/keying/combine. gpu_ = std::make_unique(); - gpu_active_ = gpu_->Available(); + gpu_active_ = gpu_->Available() && (std::getenv("JFJOCH_RSM_NO_GPU") == nullptr); if (gpu_active_) { const int n = static_cast(partials.size()); - std::vector I(n), sigma(n), rlp(n), part(n), zeta(n), corr(n); + std::vector I(n), sigma(n), rlp(n), part(n), zeta(n), corr(n), bkg(n), img(n), dd(n); std::vector onice(n); std::vector frm(n); for (int i = 0; i < n; ++i) { const auto &o = partials[i]; I[i] = o.I; sigma[i] = o.sigma; rlp[i] = o.rlp; part[i] = o.partiality; zeta[i] = o.zeta; onice[i] = o.on_ice; frm[i] = o.frame; corr[i] = o.corr; + bkg[i] = o.bkg; img[i] = o.image_number; dd[i] = o.d; } gpu_->SetPartials(n, n_frames, I.data(), sigma.data(), rlp.data(), part.data(), zeta.data(), onice.data(), frm.data(), corr.data(), frame_start.data(), frame_count.data()); - logger.Info("RotationScaleMerge: GPU partial-scaling active"); + gpu_->SetCombineInputs(bkg.data(), img.data(), dd.data()); + gpu_->SetRawRuns(static_cast(rawrun_start.size()), static_cast(perm.size()), perm.data(), + rawrun_start.data(), rawrun_count.data(), + rawrun_h.data(), rawrun_k.data(), rawrun_l.data()); + gpu_combine_ = std::getenv("JFJOCH_RSM_GPU_COMBINE") != nullptr; + logger.Info("RotationScaleMerge: GPU partial-scaling active{}", + gpu_combine_ ? " (+ GPU combine)" : ""); } #endif } @@ -643,6 +651,11 @@ void RotationScaleMerge::Combine() { fulls.insert(fulls.end(), part[t].begin(), part[t].end()); } + SortFullsByFrame(); + logger.Info("3D combine: {} fulls from {} partials", fulls.size(), n_used); +} + +void RotationScaleMerge::SortFullsByFrame() { // Sort the fulls by their (peak) frame and build per-frame CSR ranges, so the scale-fulls step can // fit a per-frame G by slicing contiguous ranges (the same layout the partials use). std::sort(fulls.begin(), fulls.end(), @@ -657,7 +670,6 @@ void RotationScaleMerge::Combine() { fulls_frame_count[f] = j - i; i = j; } - logger.Info("3D combine: {} fulls from {} partials", fulls.size(), n_used); } void RotationScaleMerge::FinalizePerFrameScale(int n_groups, const std::vector &partial_group_mean, @@ -1085,7 +1097,37 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, FinalizePerFrameScale(n_groups, partial_mean, partial_scaled); // --- 3. 3D combine of per-frame partials into fulls (fulls inherit their ASU group here). --- - Combine(); + bool combined_on_gpu = false; +#ifdef JFJOCH_USE_CUDA + // The GPU combine mirrors Combine() exactly but keeps the fulls on the device; the diagnostic dump + // (serial, one writer) has no GPU equivalent, so fall back to the CPU path when it is requested. + if (gpu_active_ && gpu_combine_ && observation_dump_path.empty()) { + std::vector corr(partials.size()); // refresh the smoothed corr on the device + for (size_t i = 0; i < partials.size(); ++i) corr[i] = partials[i].corr; + gpu_->SetCorr(corr.data()); + const int nf = gpu_->Combine(rawrun_group.data(), min_partiality, capture_uncertainty_coeff); + fulls.assign(nf, Obs{}); + g_full.assign(n_frames, 1.0); + std::vector fh(nf), fk(nf), fl(nf), fframe(nf), fgroup(nf); + std::vector fI(nf), fsig(nf), fd(nf), fimg(nf); + std::vector fon(nf); + gpu_->GetFulls(fh.data(), fk.data(), fl.data(), fI.data(), fsig.data(), fd.data(), + fimg.data(), fframe.data(), fon.data(), fgroup.data()); + for (int i = 0; i < nf; ++i) { + Obs &o = fulls[i]; + o.h = fh[i]; o.k = fk[i]; o.l = fl[i]; + o.I = fI[i]; o.sigma = fsig[i]; o.d = fd[i]; + o.rlp = 1.0f; o.partiality = 1.0f; o.corr = 1.0f; + o.image_number = fimg[i]; o.frame = fframe[i]; + o.on_ice = fon[i]; o.group = fgroup[i]; + } + SortFullsByFrame(); + logger.Info("3D combine (GPU): {} fulls", nf); + combined_on_gpu = true; + } +#endif + if (!combined_on_gpu) + Combine(); lap("combine"); // --- 4. Scale the fulls (XDS order, Unity model). --- diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index d858a51a..aba7911f 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -133,6 +133,11 @@ private: // inactive when no GPU; the CPU loops are the fallback. Built in Ingest. std::unique_ptr gpu_; bool gpu_active_ = false; + // Phase-2 GPU combine (partials->fulls on the device). Opt-in for now (JFJOCH_RSM_GPU_COMBINE): + // it is validated bit-for-bit against the CPU combine and deterministic, but until the fulls stay + // resident for the scale-fulls + merge stages the fulls round-trip makes it no faster than the CPU + // path, so it is not on by default. See RotationScaleMergeGPU::Combine. + bool gpu_combine_ = false; #endif // --- helpers (each a flat pass; see the .cpp) --- @@ -165,7 +170,11 @@ private: // that prediction applied. SG-independent, so done once in Ingest. void SmoothMosaicityAndPartiality(); - void Combine(); // partials -> fulls + void Combine(); // partials -> fulls (CPU) + + // Sort `fulls` by peak frame and (re)build fulls_frame_start/count (the per-frame CSR the scale-fulls + // step slices). Shared by the CPU Combine tail and the GPU combine path. + void SortFullsByFrame(); // Per-frame CC vs the partial merge reference, then write G/CC/mosaicity back onto the partials // (once, at the end of partial scaling) so the offline per-image scaling table is still exported. diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index b0e0eb9b..3419913d 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -174,6 +174,168 @@ namespace { } } + // std::max / std::min return (a0, sigma>0. + __device__ __forceinline__ bool CombineUsable(int i, const float *I, const float *sigma, + const float *corr) { + const float c = corr[i]; + if (!(c > 0.0f) || !isfinite(c)) return false; + return isfinite(I[i]) && isfinite(sigma[i]) && sigma[i] > 0.0f; + } + + // All device pointers + scalars the combine kernels need, passed by value. + struct CombineParams { + int n_runs; + double min_partiality, capture_uncertainty_coeff; + const float *I, *sigma, *corr, *partiality, *bkg, *image_number, *d; + const int32_t *frame; + const uint8_t *on_ice; + const int32_t *perm, *rr_start, *rr_count, *rr_h, *rr_k, *rr_l, *rr_group; + int32_t *rr_nevents, *rr_nusable; // count pass outputs + const int32_t *rr_offset; // emit pass: per-run base offset into the fulls arrays + int32_t *f_h, *f_k, *f_l, *f_frame, *f_group; + float *f_I, *f_sigma, *f_d, *f_img; + uint8_t *f_on_ice; + }; + + // One thread's raw-hkl run: split its usable partials (already in image_number order within the run) + // into rocking events, and for each event pool background, seed F and run the 3-iter de-biased Poisson + // reweight - the exact objective of RotationScaleMerge::Combine::process_rawrun. When Emit, write the + // resulting full at rr_offset[r] + (event index); otherwise just count the emitted events. Both modes + // run the identical accept test, so the count pass predicts the emit pass exactly. + template + __device__ void CombineRawRun(int r, const CombineParams &p) { + const int lo = p.rr_start[r]; + const int hi = lo + p.rr_count[r]; + const int group = p.rr_group[r]; + + int n_emit = 0; + int cursor = lo; + while (cursor < hi) { + while (cursor < hi && !CombineUsable(p.perm[cursor], p.I, p.sigma, p.corr)) ++cursor; + if (cursor >= hi) break; + const int ev_start = cursor; // first usable position of the event + int ev_end = cursor; // last usable position (inclusive), extended below + float last_img = p.image_number[p.perm[cursor]]; + int probe = cursor + 1; + while (probe < hi) { + while (probe < hi && !CombineUsable(p.perm[probe], p.I, p.sigma, p.corr)) ++probe; + if (probe >= hi) break; + const float img = p.image_number[p.perm[probe]]; + if (img - last_img > COMBINE_MAX_FRAME_GAP) break; + last_img = img; + ev_end = probe; + ++probe; + } + cursor = ev_end + 1; + + // Pass A: pooled background = mean of the event members' finite backgrounds. + double pooled_bkg = 0.0; + int n_pool = 0; + for (int m = ev_start; m <= ev_end; ++m) { + const int i = p.perm[m]; + if (!CombineUsable(i, p.I, p.sigma, p.corr)) continue; + const float b = p.bkg[i]; + if (isfinite(b)) { pooled_bkg += b; ++n_pool; } + } + pooled_bkg = n_pool > 0 ? pooled_bkg / n_pool : 0.0; + + auto pooled_I = [&](int i) -> double { + const double n_bkg = Dmax(0.0, double(p.sigma[i]) * p.sigma[i] - p.I[i]) + / Fmax(p.bkg[i], 1.0f); + return double(p.I[i]) + n_bkg * (double(p.bkg[i]) - pooled_bkg); + }; + + // Pass B: seed F (inverse-variance mean of pooled_I*corr), plus peak / d / on_ice / partiality. + double sum_w = 0.0, sum_wI = 0.0, sum_partiality = 0.0; + float d = NAN; + const int first = p.perm[ev_start]; + int peak_outcome = p.frame[first]; + float peak_frame = p.image_number[first]; + float peak_partiality = -1.0f; + const bool on_ice = p.on_ice[first]; + for (int m = ev_start; m <= ev_end; ++m) { + const int i = p.perm[m]; + if (!CombineUsable(i, p.I, p.sigma, p.corr)) continue; + const double sigma_corr = double(p.sigma[i]) * p.corr[i]; + const double w = 1.0 / (sigma_corr * sigma_corr); + sum_w += w; + sum_wI += w * pooled_I(i) * p.corr[i]; + sum_partiality += p.partiality[i]; + if (p.partiality[i] > peak_partiality) { + peak_partiality = p.partiality[i]; + peak_outcome = p.frame[i]; + peak_frame = p.image_number[i]; + } + if (!isfinite(d) && isfinite(p.d[i]) && p.d[i] > 0.0f) d = p.d[i]; + } + double F = sum_wI / sum_w; + + // Pass C: 3 de-biased Poisson reweights (variance = bkg part + corr*max(0,F)). + for (int iter = 0; iter < 3; ++iter) { + sum_w = 0.0; sum_wI = 0.0; + for (int m = ev_start; m <= ev_end; ++m) { + const int i = p.perm[m]; + if (!CombineUsable(i, p.I, p.sigma, p.corr)) continue; + const double corr = p.corr[i]; + const double I_corr = pooled_I(i) * corr; + const double sigma_corr = double(p.sigma[i]) * corr; + const double bkg_var = sigma_corr * sigma_corr - corr * I_corr; + double var = Dmax(0.0, bkg_var) + corr * Dmax(0.0, F); + if (!(var > 0.0)) var = sigma_corr * sigma_corr; + const double w = 1.0 / var; + sum_w += w; + sum_wI += w * I_corr; + } + F = sum_wI / sum_w; + } + + if (sum_w <= 0.0 || sum_partiality < p.min_partiality) + continue; + + double sigma_full = 1.0 / sqrt(sum_w); + if (p.capture_uncertainty_coeff > 0.0) { + const double frac = Dmin(1.0, sum_partiality); + const double extra = p.capture_uncertainty_coeff * (1.0 - frac) * Dmax(0.0, F); + sigma_full = sqrt(sigma_full * sigma_full + extra * extra); + } + + if (Emit) { + const int o = p.rr_offset[r] + n_emit; + p.f_h[o] = p.rr_h[r]; p.f_k[o] = p.rr_k[r]; p.f_l[o] = p.rr_l[r]; + p.f_I[o] = float(F); + p.f_sigma[o] = float(sigma_full); + p.f_d[o] = d; + p.f_img[o] = peak_frame; + p.f_frame[o] = peak_outcome; + p.f_on_ice[o] = on_ice ? 1 : 0; + p.f_group[o] = group; + } + ++n_emit; + } + + if (!Emit) { + p.rr_nevents[r] = n_emit; + int n_usable = 0; + for (int m = lo; m < hi; ++m) + if (CombineUsable(p.perm[m], p.I, p.sigma, p.corr)) ++n_usable; + p.rr_nusable[r] = n_usable; + } + } + + template + __global__ void CombineKernel(CombineParams p) { + for (int r = blockIdx.x * blockDim.x + threadIdx.x; r < p.n_runs; r += gridDim.x * blockDim.x) + CombineRawRun(r, p); + } + void CudaCheck(cudaError_t e, const char *what) { if (e != cudaSuccess) throw JFJochException(JFJochExceptionCategory::GPUCUDAError, @@ -205,6 +367,17 @@ struct RotationScaleMergeGPU::Impl { CudaDevicePtr scaled; CudaDevicePtr sco_coeff; CudaDevicePtr sco_ok; + + // combine: extra per-obs inputs + the one-time raw-hkl run layout + CudaDevicePtr bkg, image_number, d_obs; + int n_runs = 0, n_perm = 0; + CudaDevicePtr perm, rr_start, rr_count, rr_h, rr_k, rr_l, rr_group; + CudaDevicePtr rr_nevents, rr_nusable, rr_offset; + // combine: resident fulls SoA (rebuilt each Combine) + int n_fulls = 0; + CudaDevicePtr f_h, f_k, f_l, f_frame, f_group; + CudaDevicePtr f_I, f_sigma, f_d, f_img; + CudaDevicePtr f_on_ice; }; RotationScaleMergeGPU::RotationScaleMergeGPU() : impl_(std::make_unique()) { @@ -287,3 +460,100 @@ void RotationScaleMergeGPU::GetG(double *g_out, uint8_t *scaled_out) const { CudaCheck(cudaMemcpy(scaled_out, impl_->scaled.get(), size_t(impl_->n_frames) * sizeof(uint8_t), cudaMemcpyDeviceToHost), "download scaled"); } + +void RotationScaleMergeGPU::SetCombineInputs(const float *bkg, const float *image_number, const float *d) { + auto &dd = *impl_; + Upload(dd.bkg, bkg, dd.n_obs); + Upload(dd.image_number, image_number, dd.n_obs); + Upload(dd.d_obs, d, dd.n_obs); +} + +void RotationScaleMergeGPU::SetRawRuns(int n_runs, int n_perm, const int32_t *perm, + const int32_t *rr_start, const int32_t *rr_count, + const int32_t *rr_h, const int32_t *rr_k, const int32_t *rr_l) { + auto &d = *impl_; + d.n_runs = n_runs; + d.n_perm = n_perm; + Upload(d.perm, perm, n_perm); + Upload(d.rr_start, rr_start, n_runs); + Upload(d.rr_count, rr_count, n_runs); + Upload(d.rr_h, rr_h, n_runs); + Upload(d.rr_k, rr_k, n_runs); + Upload(d.rr_l, rr_l, n_runs); + d.rr_group = CudaDevicePtr(std::max(1, n_runs)); + d.rr_nevents = CudaDevicePtr(std::max(1, n_runs)); + d.rr_nusable = CudaDevicePtr(std::max(1, n_runs)); + d.rr_offset = CudaDevicePtr(std::max(1, n_runs)); +} + +int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_partiality, + double capture_uncertainty_coeff) { + auto &d = *impl_; + CudaCheck(cudaMemcpy(d.rr_group.get(), rawrun_group, size_t(d.n_runs) * sizeof(int32_t), + cudaMemcpyHostToDevice), "upload rr_group"); + + CombineParams p{}; + p.n_runs = d.n_runs; + p.min_partiality = min_partiality; + p.capture_uncertainty_coeff = capture_uncertainty_coeff; + p.I = d.I.get(); p.sigma = d.sigma.get(); p.corr = d.corr.get(); p.partiality = d.partiality.get(); + p.bkg = d.bkg.get(); p.image_number = d.image_number.get(); p.d = d.d_obs.get(); + p.frame = d.frame.get(); p.on_ice = d.on_ice.get(); + p.perm = d.perm.get(); p.rr_start = d.rr_start.get(); p.rr_count = d.rr_count.get(); + p.rr_h = d.rr_h.get(); p.rr_k = d.rr_k.get(); p.rr_l = d.rr_l.get(); p.rr_group = d.rr_group.get(); + p.rr_nevents = d.rr_nevents.get(); p.rr_nusable = d.rr_nusable.get(); + + const int blocks = std::min(65535, (d.n_runs + BLK - 1) / BLK); + + // Count pass: how many fulls each run emits. + CombineKernel<<>>(p); + CudaCheck(cudaGetLastError(), "combine count launch"); + + // Exclusive prefix sum on the host (deterministic) -> per-run output offset + total fulls. + std::vector nevents(d.n_runs); + CudaCheck(cudaMemcpy(nevents.data(), d.rr_nevents.get(), size_t(d.n_runs) * sizeof(int32_t), + cudaMemcpyDeviceToHost), "download nevents"); + std::vector offset(d.n_runs); + int64_t acc = 0; + for (int r = 0; r < d.n_runs; ++r) { offset[r] = static_cast(acc); acc += nevents[r]; } + d.n_fulls = static_cast(acc); + + // Allocate the fulls SoA and emit. + const int nf = std::max(1, d.n_fulls); + d.f_h = CudaDevicePtr(nf); d.f_k = CudaDevicePtr(nf); d.f_l = CudaDevicePtr(nf); + d.f_frame = CudaDevicePtr(nf); d.f_group = CudaDevicePtr(nf); + d.f_I = CudaDevicePtr(nf); d.f_sigma = CudaDevicePtr(nf); + d.f_d = CudaDevicePtr(nf); d.f_img = CudaDevicePtr(nf); + d.f_on_ice = CudaDevicePtr(nf); + CudaCheck(cudaMemcpy(d.rr_offset.get(), offset.data(), size_t(d.n_runs) * sizeof(int32_t), + cudaMemcpyHostToDevice), "upload offset"); + + p.rr_offset = d.rr_offset.get(); + p.f_h = d.f_h.get(); p.f_k = d.f_k.get(); p.f_l = d.f_l.get(); + p.f_frame = d.f_frame.get(); p.f_group = d.f_group.get(); + p.f_I = d.f_I.get(); p.f_sigma = d.f_sigma.get(); p.f_d = d.f_d.get(); p.f_img = d.f_img.get(); + p.f_on_ice = d.f_on_ice.get(); + if (d.n_fulls > 0) { + CombineKernel<<>>(p); + CudaCheck(cudaGetLastError(), "combine emit launch"); + } + CudaCheck(cudaDeviceSynchronize(), "combine sync"); + return d.n_fulls; +} + +void RotationScaleMergeGPU::GetFulls(int32_t *h, int32_t *k, int32_t *l, float *I, float *sigma, float *d, + float *image_number, int32_t *frame, uint8_t *on_ice, + int32_t *group) const { + const auto &dd = *impl_; + const size_t n = static_cast(dd.n_fulls); + if (n == 0) return; + auto dl = [&](void *dst, const void *src, size_t bytes) { + CudaCheck(cudaMemcpy(dst, src, bytes, cudaMemcpyDeviceToHost), "download fulls"); + }; + dl(h, dd.f_h.get(), n * sizeof(int32_t)); dl(k, dd.f_k.get(), n * sizeof(int32_t)); + dl(l, dd.f_l.get(), n * sizeof(int32_t)); dl(frame, dd.f_frame.get(), n * sizeof(int32_t)); + dl(group, dd.f_group.get(), n * sizeof(int32_t)); + dl(I, dd.f_I.get(), n * sizeof(float)); dl(sigma, dd.f_sigma.get(), n * sizeof(float)); + dl(d, dd.f_d.get(), n * sizeof(float)); dl(image_number, dd.f_img.get(), n * sizeof(float)); + dl(on_ice, dd.f_on_ice.get(), n * sizeof(uint8_t)); +} diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index cb01e49f..7963b8b9 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -54,6 +54,32 @@ public: void GetCorr(float *corr_out) const; void GetG(double *g_out, uint8_t *scaled_out) const; + // --- 3D combine (partials -> fulls), all on the device --- + + // The per-obs fields the combine needs on top of the scaling inputs (image-local bkg, fractional + // frame position for event contiguity, resolution). Uploaded once, alongside SetPartials. + void SetCombineInputs(const float *bkg, const float *image_number, const float *d); + + // The one-time raw-hkl run layout (space-group-independent): the (raw h,k,l, image_number)-sorted + // permutation of the obs, split into contiguous per-raw-hkl runs. Uploaded once in Ingest. + void SetRawRuns(int n_runs, int n_perm, const int32_t *perm, + const int32_t *rawrun_start, const int32_t *rawrun_count, + const int32_t *rawrun_h, const int32_t *rawrun_k, const int32_t *rawrun_l); + + // Combine the resident partials (reading the current resident corr) into fulls on the device, + // mirroring RotationScaleMerge::Combine: one thread per raw-hkl run splits its usable partials into + // rocking events (frame gap <= 2), pools background, seeds F, does 3 de-biased Poisson reweights and + // adds the capture-uncertainty term. rawrun_group (length n_runs) is the current space group's ASU + // id per raw hkl (it becomes the full's group). Deterministic: fulls are emitted in raw-run-major, + // event order (a count pass -> host prefix sum -> emit-at-offset), matching the CPU path. Returns the + // number of fulls (call GetFulls with buffers of that length). + int Combine(const int32_t *rawrun_group, double min_partiality, double capture_uncertainty_coeff); + + // Download the combined fulls SoA (length = Combine()'s return). corr/partiality/rlp are 1 by + // construction and not returned; the caller sets them. + void GetFulls(int32_t *h, int32_t *k, int32_t *l, float *I, float *sigma, float *d, + float *image_number, int32_t *frame, uint8_t *on_ice, int32_t *group) const; + private: struct Impl; std::unique_ptr impl_;