mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases 67 Wiki Activity

RotationScaleMerge: GPU 3D-combine (CUDA port, phase 2 step 1)

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-03 07:27:33 +02:00
co-authored by Claude Opus 4.8
parent dd25de461d
commit 2c928b27cd
4 changed files with 353 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+47 -5
View File
@@ -7,6 +7,7 @@
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <fstream>
#include <future>
#include <limits>
@@ -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<RotationScaleMergeGPU>();
gpu_active_ = gpu_->Available();
gpu_active_ = gpu_->Available() && (std::getenv("JFJOCH_RSM_NO_GPU") == nullptr);
if (gpu_active_) {
const int n = static_cast<int>(partials.size());
std::vector<float> I(n), sigma(n), rlp(n), part(n), zeta(n), corr(n);
std::vector<float> I(n), sigma(n), rlp(n), part(n), zeta(n), corr(n), bkg(n), img(n), dd(n);
std::vector<uint8_t> onice(n);
std::vector<int32_t> 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<int>(rawrun_start.size()), static_cast<int>(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<double> &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<float> 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<int32_t> fh(nf), fk(nf), fl(nf), fframe(nf), fgroup(nf);
std::vector<float> fI(nf), fsig(nf), fd(nf), fimg(nf);
std::vector<uint8_t> 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). ---
image_analysis/scale_merge/RotationScaleMerge.h
+10 -1
View File
@@ -133,6 +133,11 @@ private:
// inactive when no GPU; the CPU loops are the fallback. Built in Ingest.
std::unique_ptr<RotationScaleMergeGPU> 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.
image_analysis/scale_merge/RotationScaleMergeGPU.cu
+270
View File
@@ -174,6 +174,168 @@ namespace {
}
}
// std::max / std::min return (a<b)?b:a and (b<a)?b:a - reproduce that exactly (NOT fmax/fmin, which
// differ on NaN) so the combine matches the CPU path bit-for-bit on the same inputs.
__device__ __forceinline__ double Dmax(double a, double b) { return (a < b) ? b : a; }
__device__ __forceinline__ double Dmin(double a, double b) { return (b < a) ? b : a; }
__device__ __forceinline__ float Fmax(float a, float b) { return (a < b) ? b : a; }
constexpr float COMBINE_MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge::MAX_FRAME_GAP
// A partial is usable for the combine iff its corr and (I, sigma) are finite with corr>0, 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 <bool Emit>
__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 <bool Emit>
__global__ void CombineKernel(CombineParams p) {
for (int r = blockIdx.x * blockDim.x + threadIdx.x; r < p.n_runs; r += gridDim.x * blockDim.x)
CombineRawRun<Emit>(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<uint8_t> scaled;
CudaDevicePtr<float> sco_coeff;
CudaDevicePtr<uint8_t> sco_ok;
// combine: extra per-obs inputs + the one-time raw-hkl run layout
CudaDevicePtr<float> bkg, image_number, d_obs;
int n_runs = 0, n_perm = 0;
CudaDevicePtr<int32_t> perm, rr_start, rr_count, rr_h, rr_k, rr_l, rr_group;
CudaDevicePtr<int32_t> rr_nevents, rr_nusable, rr_offset;
// combine: resident fulls SoA (rebuilt each Combine)
int n_fulls = 0;
CudaDevicePtr<int32_t> f_h, f_k, f_l, f_frame, f_group;
CudaDevicePtr<float> f_I, f_sigma, f_d, f_img;
CudaDevicePtr<uint8_t> f_on_ice;
};
RotationScaleMergeGPU::RotationScaleMergeGPU() : impl_(std::make_unique<Impl>()) {
@@ -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<int32_t>(std::max(1, n_runs));
d.rr_nevents = CudaDevicePtr<int32_t>(std::max(1, n_runs));
d.rr_nusable = CudaDevicePtr<int32_t>(std::max(1, n_runs));
d.rr_offset = CudaDevicePtr<int32_t>(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<false><<<blocks, BLK>>>(p);
CudaCheck(cudaGetLastError(), "combine count launch");
// Exclusive prefix sum on the host (deterministic) -> per-run output offset + total fulls.
std::vector<int32_t> 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<int32_t> offset(d.n_runs);
int64_t acc = 0;
for (int r = 0; r < d.n_runs; ++r) { offset[r] = static_cast<int32_t>(acc); acc += nevents[r]; }
d.n_fulls = static_cast<int>(acc);
// Allocate the fulls SoA and emit.
const int nf = std::max(1, d.n_fulls);
d.f_h = CudaDevicePtr<int32_t>(nf); d.f_k = CudaDevicePtr<int32_t>(nf); d.f_l = CudaDevicePtr<int32_t>(nf);
d.f_frame = CudaDevicePtr<int32_t>(nf); d.f_group = CudaDevicePtr<int32_t>(nf);
d.f_I = CudaDevicePtr<float>(nf); d.f_sigma = CudaDevicePtr<float>(nf);
d.f_d = CudaDevicePtr<float>(nf); d.f_img = CudaDevicePtr<float>(nf);
d.f_on_ice = CudaDevicePtr<uint8_t>(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<true><<<blocks, BLK>>>(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<size_t>(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));
}
image_analysis/scale_merge/RotationScaleMergeGPU.h
+26
View File
@@ -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> impl_;