From f46dc303a3ecd235d00fd824c1835beb47413fa7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 15 Oct 2026 11:13:41 +0200 Subject: [PATCH] Rotation scaling on the GPU: two-pass ingest, decided before the card is touched, output unchanged - The resident ingest reads the reflections twice: first only what the sort and the cut need (frame position, resolution, background, raw-hkl key), then - once the GPU/CPU decision has been made on the exact counts - the observations the cut kept, written straight to their place. The observation arrays and the sort buffers are never on the card together, and no observation field is cut (held twice) any more. - A set too large even for that sort scales on the CPU before anything is allocated on the card (ScaleOnCPUWhereSortTooLarge, same 80%-of-total rule). - The predicted detector position never goes to the card: it is taken on the host, and the combine emits the index of each full's peak partial instead of its px/py (8 B less per observation and 4 per full). - A full carries its raw-hkl run index instead of h, k, l (8 B less per full); the constant-1 partiality array of the fulls is gone. - The fulls' CSR sorts and the error-model rank sort run on cub::DoubleBuffer. Measured with per-array device accounting on an RTX 5080, largest engine of each set, before (c130c7923) -> after, MB: set obs/hkl ingest peak (B/refl) post-decision peak 4nwv 32.7 M/15.4 M 3277 (100) -> 2336 (72) 4265 -> 3625 8qaw 122.7 M/ 8.6 M 10756 (88) -> 7421 (61) 9327 -> 8156 9v16 54.7 M/ 5.7 M 4838 (89) -> 3419 (63) 4283 -> 3774 7qij 78.1 M/ 2.2 M 6768 (87) -> 4707 (60) 5242 -> 4707 8sa8 50.0 M/ 5.0 M 4107 (82) -> 3134 (63) 4205 -> 3603 The sort before the decision now holds 32 B per reflection (46 on a P1-like set with many runs); the engine's overall peak is after the decision on all of them. GPUScalingPeakBytes refitted on the same fourteen engines: 63 + 115 -> 56 B/obs + 99 B/hkl, within 9%. DeviceBytes as an upper bound: 66/94/52 -> 58/30/86 per obs/event/hkl, ingest 89/56 -> 60/42 per reflection/run. p.hkl byte-identical to the rc175 baseline on 15 battery sets (4nwv, 8qaw, 9v16, 7qij, 8sa8, 11if, 5reo, 6oel, 9jzo, 9min, metformin, lyso_x06da_ref, thau_bl1a_3p8keV, insu_I_x06da_5keV_2, myob_x06da); compare shows no change. [RotationScale] (RotationIngestGPU now also checks the host positions) and ErrorModel_DeviceFitIsTheHostFit pass. Co-Authored-By: Claude Opus 5.5 --- image_analysis/scale_merge/ErrorModelGPU.cu | 21 +- .../scale_merge/RotationScaleMerge.cpp | 105 +++-- .../scale_merge/RotationScaleMerge.h | 20 +- .../scale_merge/RotationScaleMergeGPU.cu | 368 +++++++++++------- .../scale_merge/RotationScaleMergeGPU.h | 42 +- rugnux/RugnuxScaleMerge.cpp | 2 +- tests/RotationIngestGPUTest.cpp | 10 +- 7 files changed, 346 insertions(+), 222 deletions(-) diff --git a/image_analysis/scale_merge/ErrorModelGPU.cu b/image_analysis/scale_merge/ErrorModelGPU.cu index ef1302b34..8d56d5631 100644 --- a/image_analysis/scale_merge/ErrorModelGPU.cu +++ b/image_analysis/scale_merge/ErrorModelGPU.cu @@ -89,8 +89,8 @@ ErrorModelFit FitErrorModelGPU(const std::vector &pool) { return fit; const int n = static_cast(pool.size()); - // The rank sort is the peak: the pool, two key and two index arrays and the sort scratch. - constexpr size_t BYTES_PER_SAMPLE = 32 + 2 * 8 + 2 * 4 + 16; + // The rank sort is the peak: the pool, two key and two index arrays. + constexpr size_t BYTES_PER_SAMPLE = 32 + 2 * 8 + 2 * 4; size_t free_bytes = 0; if (!cuda_memory_available(BYTES_PER_SAMPLE * pool.size(), free_bytes)) throw JFJochException(JFJochExceptionCategory::GPUCUDAError, @@ -104,22 +104,25 @@ ErrorModelFit FitErrorModelGPU(const std::vector &pool) { stream)); // The samples in the host's rank order - (I2/s2, s2, I2, dev2), a total order - by four stable - // passes, least significant key first. The equal-count bins are then consecutive stretches. + // passes, least significant key first. The equal-count bins are then consecutive stretches. Each + // pass goes back and forth between the two key and the two index arrays, so that cub holds no copy + // of them of its own; the same stable sort, so the same order. CudaDevicePtr idx(n, ALLOC), idx_next(n, ALLOC); CudaDevicePtr key(n, ALLOC), key_sorted(n, ALLOC); IotaKernel<<>>(n, idx); cuda_err(cudaGetLastError()); + cub::DoubleBuffer keys(key.get(), key_sorted.get()); + cub::DoubleBuffer order(idx.get(), idx_next.get()); size_t bytes = 0; - cuda_err(cub::DeviceRadixSort::SortPairs(nullptr, bytes, key.get(), key_sorted.get(), idx.get(), idx_next.get(), - n, 0, 64, stream)); + cuda_err(cub::DeviceRadixSort::SortPairs(nullptr, bytes, keys, order, n, 0, 64, stream)); CudaDevicePtr scratch(std::max(bytes, 1), ALLOC); for (int field = 0; field < 4; ++field) { - KeyKernel<<>>(n, d_pool, idx, field, key); + KeyKernel<<>>(n, d_pool, order.Current(), field, keys.Current()); cuda_err(cudaGetLastError()); - cuda_err(cub::DeviceRadixSort::SortPairs(scratch.get(), bytes, key.get(), key_sorted.get(), idx.get(), - idx_next.get(), n, 0, 64, stream)); - std::swap(idx, idx_next); + cuda_err(cub::DeviceRadixSort::SortPairs(scratch.get(), bytes, keys, order, n, 0, 64, stream)); } + if (keys.Current() != key_sorted.get()) std::swap(key, key_sorted); + if (order.Current() != idx.get()) std::swap(idx, idx_next); scratch = CudaDevicePtr(); idx_next = CudaDevicePtr(); key = CudaDevicePtr(); // After the last pass key_sorted holds I2/s2 in rank order. CudaDevicePtr s2(n, ALLOC), I2(n, ALLOC), dev2(n, ALLOC); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7595dcca8..7ab124f79 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -608,18 +608,23 @@ void RotationScaleMerge::IngestHost() { } #ifdef JFJOCH_USE_CUDA -// The resident path's ingest. The reflections go to the card once, frame by frame, and the card measures -// the per-frame background, sorts, cuts the raw-hkl runs and drops what the resolution range excludes -// (RotationScaleMergeGPU::IngestUpload .. IngestSort) - the observations, the runs and the order of each -// that BuildInRangeObservations lays down on the host, to the bit, and left where the device pipeline -// reads them. The geometry smoothing then runs here on what the card hands back. +// The resident path's ingest. The reflections go to the card frame by frame, and the card measures the +// per-frame background, sorts, cuts the raw-hkl runs and drops what the resolution range excludes +// (RotationScaleMergeGPU::IngestUpload .. IngestSort); once the scaling is known to fit, they go to the +// card a second time for the observations the cut kept (IngestObservations) - the observations, the runs +// and the order of each that BuildInRangeObservations lays down on the host, to the bit, and left where +// the device pipeline reads them. The geometry smoothing then runs here on what the card hands back. void RotationScaleMerge::IngestOnDevice() { std::vector refl(n_frames); std::vector count(n_frames); + size_t n_reflections = 0; for (int o = 0; o < n_frames; ++o) { refl[o] = partials_out[o].reflections.data(); count[o] = static_cast(partials_out[o].reflections.size()); + n_reflections += count[o]; } + if (ScaleOnCPUWhereSortTooLarge(n_reflections)) + return; TimingMark upload_mark(*logger, "rsm ingest upload"); gpu_->IngestUpload(n_frames, refl.data(), count.data(), nthreads); upload_mark.Stop(); @@ -627,7 +632,7 @@ void RotationScaleMerge::IngestOnDevice() { std::vector mean_bkg(n_frames); gpu_->IngestMeanBackground(mean_bkg.data()); MeasureIncidentFlux(mean_bkg); - const int n = gpu_->IngestSort(frame_flux.data(), d_min_limit, d_max_limit); + const int n = gpu_->IngestSort(d_min_limit, d_max_limit); n_partials_obs = n; RotationScaleMergeGPU::IngestLayout layout; gpu_->IngestGetLayout(layout); @@ -644,9 +649,13 @@ void RotationScaleMerge::IngestOnDevice() { rawrun_l = std::move(layout.run_l); rawrun_d = std::move(layout.run_d); // The scaling's footprint is known from here on: where the card cannot take it, the card's work - // ends before the smoothing. + // ends before a single observation is on it. if (ScaleOnCPUWhereTooLarge()) return; + auto pos = std::make_shared(); + pos->px.resize(n); pos->py.resize(n); + gpu_->IngestObservations(refl.data(), count.data(), nthreads, frame_flux.data(), pos->px.data(), pos->py.data()); + obs_position = std::move(pos); const SmoothedGeometry sg = SmoothedFrameGeometry(); int64_t n_resolved = 0; @@ -703,13 +712,28 @@ void RotationScaleMerge::IngestOnDevice() { // scales on the CPU - decided here, once, on the card's total memory, so the answer depends on the card // but not on the moment. bool RotationScaleMerge::ScaleOnCPUWhereTooLarge() { + return ScaleOnCPUWhereOverCard(GPUScalingPeakBytes(), fmt::format("scaling {} observations of {} distinct hkl", + n_partials_obs, rawrun_start.size())); +} + +// Before the sort no observation is on the card yet, nor is it known how many there are: what the sort +// and the cut hold (IngestUpload .. IngestSort), measured at 32 bytes per reflection, or 20 and 56 per +// raw-hkl run where the runs are many - 48 per reflection with half as many runs as reflections, the +// most measured being 0.47 - and the two staging buffers. A set that the card cannot even sort scales +// on the CPU without touching the card. +bool RotationScaleMerge::ScaleOnCPUWhereSortTooLarge(size_t n_reflections) { + constexpr size_t SORT_BYTES_PER_REFLECTION = 48, STAGING_BYTES = size_t(128) << 20; + return ScaleOnCPUWhereOverCard(SORT_BYTES_PER_REFLECTION * n_reflections + STAGING_BYTES, + fmt::format("sorting {} reflections", n_reflections)); +} + +bool RotationScaleMerge::ScaleOnCPUWhereOverCard(size_t bytes, const std::string &what) { const size_t card = get_gpu_total_memory(gpu_device); - if (GPUScalingPeakBytes() <= CPU_SCALING_ABOVE_CARD_FRACTION * card) + if (bytes <= CPU_SCALING_ABOVE_CARD_FRACTION * card) return false; - logger->Warning("RotationScaleMerge: scaling {} observations of {} distinct hkl on the GPU needs about " - "{:.1f} GB, over {:.0f}% of the card's {:.1f} GB - scaling on the CPU instead", - n_partials_obs, rawrun_start.size(), GPUScalingPeakBytes() / 1e9, - 100.0 * CPU_SCALING_ABOVE_CARD_FRACTION, card / 1e9); + logger->Warning("RotationScaleMerge: {} on the GPU needs about {:.1f} GB, over {:.0f}% of the card's {:.1f} GB " + "- scaling on the CPU instead", what, bytes / 1e9, 100.0 * CPU_SCALING_ABOVE_CARD_FRACTION, + card / 1e9); scale_on_cpu_ = true; gpu_active_ = false; return true; @@ -725,12 +749,13 @@ bool RotationScaleMerge::OnGPU() const { } size_t RotationScaleMerge::GPUScalingPeakBytes() const { - // The device scaling's peak from this decision on, fitted to the engine's device buffers measured at - // their peak on fourteen engines of five large data sets (23-123 M observations, 1.6-15 M distinct - // hkl): within 12% on each, the furthest under a P1 merge of a large cell. The ingest's own peak, 82 - // to 100 bytes per reflection, has passed by the time this decides. Unlike DeviceBytes, an estimate - // and not an upper bound - it decides where the GPU clearly cannot take the scaling. - constexpr size_t PEAK_BYTES_PER_OBS = 63, PEAK_BYTES_PER_HKL = 115; + // The device scaling's peak from this decision on - the observations taken from the reflections + // (IngestObservations) included, the sort before it not - fitted to the engine's device buffers + // measured at their peak on fourteen engines of five large data sets (23-123 M observations, + // 1.6-15 M distinct hkl): within 9% on each, the furthest under a P1 merge of a large cell. Unlike + // DeviceBytes, an estimate and not an upper bound - it decides where the GPU clearly cannot take the + // scaling. + constexpr size_t PEAK_BYTES_PER_OBS = 56, PEAK_BYTES_PER_HKL = 99; return PEAK_BYTES_PER_OBS * static_cast(n_partials_obs) + PEAK_BYTES_PER_HKL * rawrun_start.size(); } @@ -738,22 +763,21 @@ size_t RotationScaleMerge::DeviceBytes() const { // What each costs at the peak: an observation's arrays, and its group CSR until the combine; a full's // arrays with its error-model sample, or with the surface terms; an hkl's raw run and its ASU group's // sums. Fitted as an upper bound to the engine's device buffers measured at their peak on fourteen - // engines of five large data sets (23-123 M observations): up to 18% over. The events are counted + // engines of five large data sets (23-123 M observations): up to 14% over. The events are counted // on the device path's ingest; elsewhere every observation may be an event of its own. - constexpr size_t BYTES_PER_OBS = 66, BYTES_PER_EVENT = 94, BYTES_PER_HKL = 52; + constexpr size_t BYTES_PER_OBS = 58, BYTES_PER_EVENT = 30, BYTES_PER_HKL = 86; const auto events = static_cast(resident_ingest ? rocking_events_at_ingest : n_partials_obs); const size_t bytes = BYTES_PER_OBS * static_cast(n_partials_obs) + BYTES_PER_EVENT * events + BYTES_PER_HKL * rawrun_start.size(); - // The resident ingest's own peak comes before any of that, over every reflection read and every run - // before the resolution cut, measured on the same engines. Beside the observation arrays and the - // integrated delta_phi (58 bytes), each radix sort holds two pairs of keys and values, 82 bytes per - // reflection in all; the cut (IngestSort) of the runs 74 per reflection and 56 per run, that of the - // observations after it 89 per reflection and 24 per run kept. The upload before them holds the - // observation arrays, the sort key and delta_phi (66 bytes) and the two 64 MB staging buffers. - constexpr size_t INGEST_BYTES_PER_REFLECTION = 89, INGEST_BYTES_PER_RUN = 56; - constexpr size_t INGEST_ARRAY_BYTES = 66, INGEST_STAGING_BYTES = size_t(160) << 20; - const size_t ingest = std::max(INGEST_BYTES_PER_REFLECTION * ingest_reflections + INGEST_BYTES_PER_RUN * ingest_runs, - INGEST_ARRAY_BYTES * ingest_reflections + INGEST_STAGING_BYTES); + // The resident ingest's own peak, over every reflection read and every run before the resolution + // cut, measured on the same engines: the observations taken from the reflections hold the + // observation arrays, the place the cut gave each reflection and the runs it kept, 60 to 72 bytes + // per reflection; the sort and the cut before them 32 per reflection and 56 per run. The two 64 MB + // staging buffers beside either. + constexpr size_t INGEST_BYTES_PER_REFLECTION = 60, INGEST_BYTES_PER_RUN = 42; + constexpr size_t INGEST_STAGING_BYTES = size_t(160) << 20; + const size_t ingest = INGEST_BYTES_PER_REFLECTION * ingest_reflections + INGEST_BYTES_PER_RUN * ingest_runs + + INGEST_STAGING_BYTES; return std::max(bytes + bytes / 5, ingest); } @@ -799,7 +823,7 @@ void RotationScaleMerge::IngestDevice(const RotationScaleMerge *observations_fro using F = RotationScaleMergeGPU::ObsField; constexpr int STAGE_SLICE = 8 << 20; const int slice = std::min(std::max(n, 1), STAGE_SLICE); - constexpr int NFLOAT = 12; + constexpr int NFLOAT = 10; std::unique_ptr fbuf[NFLOAT]; for (auto &b : fbuf) b.reset(new float[slice]); std::unique_ptr frm(new int32_t[slice]); @@ -812,19 +836,25 @@ void RotationScaleMerge::IngestDevice(const RotationScaleMerge *observations_fro fbuf[0][i] = o.I; fbuf[1][i] = o.sigma; fbuf[2][i] = o.prescaling_corr; fbuf[3][i] = o.partiality; fbuf[4][i] = o.zeta; fbuf[5][i] = o.corr; fbuf[6][i] = o.bkg; fbuf[7][i] = o.var_bkg; fbuf[8][i] = o.image_number; - fbuf[9][i] = o.d; fbuf[10][i] = o.px; fbuf[11][i] = o.py; + fbuf[9][i] = o.d; frm[i] = o.frame; onice[i] = o.on_ice; clip[i] = o.clipped; } }); constexpr F fields[NFLOAT] = {F::I, F::Sigma, F::PrescalingCorr, F::Partiality, F::Zeta, F::Corr0, F::Bkg, F::VarBkg, F::ImageNumber, - F::D, F::Px, F::Py}; + F::D}; for (int k = 0; k < NFLOAT; ++k) gpu_->SetObsField(fields[k], at, cnt, fbuf[k].get()); gpu_->SetObsFrame(at, cnt, frm.get()); gpu_->SetObsOnIce(at, cnt, onice.get()); gpu_->SetObsClipped(at, cnt, clip.get()); } + auto pos = std::make_shared(); + pos->px.resize(n); pos->py.resize(n); + ParallelChunks(n, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) { pos->px[i] = partials[i].px; pos->py[i] = partials[i].py; } + }); + obs_position = std::move(pos); 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()); @@ -6706,23 +6736,24 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool fulls.resize(nf); FullsStaging &st = fulls_staging; st.Resize(nf); - gpu_->GetFulls(st.h.data(), st.k.data(), st.l.data(), st.I.data(), st.sigma.data(), st.d.data(), + gpu_->GetFulls(st.run.data(), st.I.data(), st.sigma.data(), st.d.data(), st.image_number.data(), st.frame.data(), st.on_ice.data(), st.clipped.data(), st.group.data()); - gpu_->GetFullsPxPy(st.px.data(), st.py.data()); + gpu_->GetFullsPeak(st.peak.data()); + const ObsPosition &pos = *obs_position; gpu_->GetFullsVariance(st.var_bkg.data(), st.var_per_I.data(), st.capture.data()); if (scaled_fulls_on_gpu) gpu_->GetFullsCorr(st.corr.data()); gpu_->FullsDownloaded(); ParallelChunks(nf, ThreadsForWork(nf, nthreads), [&](int lo, int hi) { for (int i = lo; i < hi; ++i) { Obs &o = fulls[i]; - o.h = st.h[i]; o.k = st.k[i]; o.l = st.l[i]; + o.h = rawrun_h[st.run[i]]; o.k = rawrun_k[st.run[i]]; o.l = rawrun_l[st.run[i]]; o.I = st.I[i]; o.sigma = st.sigma[i]; o.d = st.d[i]; o.var_bkg = st.var_bkg[i]; o.var_per_I = st.var_per_I[i]; o.capture = st.capture[i]; o.prescaling_corr = 1.0f; o.partiality = 1.0f; o.corr = scaled_fulls_on_gpu ? st.corr[i] : 1.0f; o.image_number = st.image_number[i]; o.frame = st.frame[i]; - o.px = st.px[i]; o.py = st.py[i]; + o.px = pos.px[st.peak[i]]; o.py = pos.py[st.peak[i]]; o.on_ice = st.on_ice[i]; o.clipped = st.clipped[i]; o.group = st.group[i]; // A full has no rocking geometry of its own and carries no background; the combine on // the device does not produce these, and a reused buffer would keep the last chain's. diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 61c13378c..514378909 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -173,8 +173,8 @@ public: // card. Call after IngestHost(). It counts what the ingest knows - the observations, the rocking // events (the combine makes at most one full of each) and the distinct hkl (at most one ASU group // each) - at what each costs at the engine's peak (the correction-surface terms, or the merge with - // its outlier test), plus a fifth. Measured on 11 data sets of 13-123 M observations, every - // engine's peak was at most 81% of this. Where the ingest on the card peaks higher, its peak. + // its outlier test), plus a fifth. Measured on 14 engines of 23-123 M observations, every engine's + // peak was at most 88% of this. Where the ingest on the card peaks higher, its peak. // Call after IngestDevice(). [[nodiscard]] size_t DeviceBytes() const; // The device scaling's estimated peak (see the .cpp). Call after the ingest. @@ -370,6 +370,10 @@ private: bool partials_released = false; int rocking_event_frames_at_ingest = 0; int64_t rocking_events_at_ingest = 0; + // Each partial's predicted detector position, for the fulls the device path combines (GetFullsPeak): + // only the host reads it, so the card does not hold it. Shared by the engines CloneIngest makes. + struct ObsPosition { std::vector px, py; }; + std::shared_ptr obs_position; // The reflections and raw-hkl runs the resident ingest read before the resolution cut (DeviceBytes). size_t ingest_reflections = 0, ingest_runs = 0; std::vector frame_start, frame_count; // CSR ranges of `partials` per frame @@ -446,13 +450,13 @@ private: // scale->combine->merge chain runs several times per run and these are a few hundred megabytes // between them, so as locals every chain allocates, faults in and zeroes the lot again. struct FullsStaging { - std::vector h, k, l, frame, group; - std::vector I, sigma, d, image_number, corr, px, py, var_bkg, var_per_I, capture; + std::vector run, peak, frame, group; + std::vector I, sigma, d, image_number, corr, var_bkg, var_per_I, capture; std::vector on_ice, clipped; void Resize(int n) { - h.resize(n); k.resize(n); l.resize(n); frame.resize(n); group.resize(n); + run.resize(n); peak.resize(n); frame.resize(n); group.resize(n); I.resize(n); sigma.resize(n); d.resize(n); image_number.resize(n); corr.resize(n); - px.resize(n); py.resize(n); var_bkg.resize(n); var_per_I.resize(n); capture.resize(n); + var_bkg.resize(n); var_per_I.resize(n); capture.resize(n); on_ice.resize(n); clipped.resize(n); } }; @@ -605,6 +609,10 @@ private: // Where the scaling's estimated peak is over CPU_SCALING_ABOVE_CARD_FRACTION of the card: says so, // gives the device up (scale_on_cpu_) and returns true. Call once the ingest has counted. bool ScaleOnCPUWhereTooLarge(); + // The same for the ingest's sort of n_reflections, decided before anything is on the card. + bool ScaleOnCPUWhereSortTooLarge(size_t n_reflections); + // What the rule above says for `bytes` of the card, `what` naming the work in the warning. + bool ScaleOnCPUWhereOverCard(size_t bytes, const std::string &what); #endif // Compute the dense ASU-group id for the current space group by grouping the (pre-sorted) raw-hkl diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index e1c425eba..612072e94 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -294,7 +294,7 @@ namespace { float max_frame_gap; const ObsRecord *__restrict__ obs; const float *__restrict__ corr, *__restrict__ bkg, *__restrict__ var_bkg, - *__restrict__ image_number, *__restrict__ d, *__restrict__ px, *__restrict__ py; + *__restrict__ image_number, *__restrict__ d; const int32_t *__restrict__ frame; const uint8_t *__restrict__ on_ice, *__restrict__ clipped; const int32_t *__restrict__ perm, *__restrict__ rr_start, *__restrict__ rr_count, @@ -303,8 +303,8 @@ namespace { int32_t *rr_nevents; // count pass output int32_t *rr_noverloaded; // count pass output: events dropped for a saturated pixel 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, *f_px, *f_py, *f_var_bkg, *f_var_per_I, *f_capture; + int32_t *f_run, *f_peak, *f_frame, *f_group; + float *f_I, *f_sigma, *f_d, *f_img, *f_var_bkg, *f_var_per_I, *f_capture; uint8_t *f_on_ice, *f_clipped; }; @@ -363,7 +363,7 @@ namespace { const int first = p.perm[ev_start]; int peak_outcome = p.frame[first]; float peak_frame = p.image_number[first]; - float peak_px = p.px[first], peak_py = p.py[first]; + int peak_obs = first; float peak_partiality = -1.0f; const bool on_ice = p.on_ice[first]; uint8_t flags = 0; @@ -383,7 +383,7 @@ namespace { peak_partiality = p.obs[i].partiality; peak_outcome = p.frame[i]; peak_frame = p.image_number[i]; - peak_px = p.px[i]; peak_py = p.py[i]; + peak_obs = i; } if (!isfinite(d) && isfinite(p.d[i]) && p.d[i] > 0.0f) d = p.d[i]; } @@ -432,7 +432,7 @@ namespace { 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_run[o] = r; p.f_I[o] = float(F); p.f_sigma[o] = float(sigma_full); p.f_var_bkg[o] = float(var_bkg_full); @@ -440,7 +440,7 @@ namespace { p.f_capture[o] = float(capture); p.f_d[o] = d; p.f_img[o] = peak_frame; - p.f_px[o] = peak_px; p.f_py[o] = peak_py; + p.f_peak[o] = peak_obs; p.f_frame[o] = peak_outcome; p.f_on_ice[o] = on_ice ? 1 : 0; p.f_clipped[o] = clipped ? 1 : 0; @@ -672,7 +672,7 @@ namespace { double min_partiality, error_model_a, error_model_b, reject_nsigma; double reject_up, reject_down; // the multiplicative part of the outlier band (OutlierBand.h) int for_search, error_model_active, reject_outliers; - const float *I, *sigma, *corr, *partiality, *d, *reject_median, *var_bkg, *var_per_I, *capture; + const float *I, *sigma, *corr, *d, *reject_median, *var_bkg, *var_per_I, *capture; const float *reject_var_add; // per group: the shell's measured share of dI^2/4 const int32_t *group, *frame; const uint8_t *on_ice, *frame_cell_ok, *half; @@ -693,6 +693,9 @@ namespace { uint8_t *rejected_obs; // per-full: set by MergeAccum (outlier-rejected), read by MergeRmeas }; + // A full's partiality: the device merges the fulls it scaled, in the Unity model. + constexpr float FULL_PARTIALITY = 1.0f; + // A full passes the merge / error-model filter (mirrors MergeAndStats::usable_merge). __device__ __forceinline__ bool MergeUsable(int i, const MergeParams &p) { const int g = p.group[i]; @@ -701,7 +704,7 @@ namespace { const float c = p.corr[i]; if (!(c > 0.0f) || !isfinite(c)) return false; if (p.for_search && p.on_ice[i]) return false; - if (p.partiality[i] < p.min_partiality) return false; + if (FULL_PARTIALITY < p.min_partiality) return false; const float I_corr = p.I[i] * c, sigma_corr = p.sigma[i] * c; return isfinite(I_corr) && isfinite(sigma_corr) && sigma_corr > 0.0f; } @@ -713,7 +716,7 @@ namespace { if (!p.frame_cell_ok[p.frame[i]]) return false; const float c = p.corr[i]; if (!(c > 0.0f) || !isfinite(c)) return false; - if (p.partiality[i] < p.min_partiality) return false; + if (FULL_PARTIALITY < p.min_partiality) return false; const float I_corr = p.I[i] * c, sigma_corr = p.sigma[i] * c; return isfinite(I_corr) && isfinite(sigma_corr) && sigma_corr > 0.0f; } @@ -838,10 +841,12 @@ namespace { } // One thread per full: its observation key, and whether the merge can use it. - __global__ void ObservationKeyKernel(int n_obs, MergeParams p, const int32_t *h, const int32_t *k, - const int32_t *l, uint64_t *key, uint8_t *usable) { + // A full's raw hkl is its run's. + __global__ void ObservationKeyKernel(int n_obs, MergeParams p, const int32_t *run, const int32_t *h, + const int32_t *k, const int32_t *l, uint64_t *key, uint8_t *usable) { for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n_obs; i += gridDim.x * blockDim.x) { - key[i] = ObservationKeyDevice(h[i], k[i], l[i], p.frame[i]); + const int r = run[i]; + key[i] = ObservationKeyDevice(h[r], k[r], l[r], p.frame[i]); usable[i] = MergeUsable(i, p) ? 1 : 0; } } @@ -1257,19 +1262,20 @@ struct RotationScaleMergeGPU::Impl { int merge_em_active = 0; // combine: extra per-obs inputs + the one-time raw-hkl run layout - CudaDevicePtr bkg, var_bkg, image_number, d_obs, px_obs, py_obs; + CudaDevicePtr bkg, var_bkg, image_number, d_obs; int n_runs = 0, n_perm = 0; CudaDevicePtr perm, rr_start, rr_count, rr_h, rr_k, rr_l; // 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, f_px, f_py, f_var_bkg, f_var_per_I, f_capture; + // run: the raw-hkl run a full is of; peak: the partial its position is taken from + CudaDevicePtr f_run, f_peak, f_frame, f_group; + CudaDevicePtr f_I, f_sigma, f_d, f_img, f_var_bkg, f_var_per_I, f_capture; CudaDevicePtr f_on_ice, f_clipped; // scale-fulls (Unity model, kept resident): all-ones partiality/zeta and records with partiality and // prescaling_corr 1, so the shared scaling kernels yield coeff=mean, plus the working corr and the // fulls frame/group CSRs (built on the host from the small f_frame/f_group key arrays, over the // emit-ordered fulls). - CudaDevicePtr f_corr, f_partiality, f_zeta; + CudaDevicePtr f_corr, f_zeta; CudaDevicePtr f_obs; CudaDevicePtr f_frame_perm, f_frame_start, f_frame_count; CudaDevicePtr f_gperm, f_gstart, f_gcount; @@ -1285,12 +1291,22 @@ struct RotationScaleMergeGPU::Impl { CudaDevicePtr s_perm[3], s_seg_start[3]; CudaDevicePtr s_A, s_sw, s_swI, s_tcross, s_tref2, s_cross, s_ref2; CudaDevicePtr s_w_Is, s_w_Iref, s_Iref; // per term of a subset, in segment order - // the resident ingest (IngestUpload .. IngestFinish): each frame's first reflection - // (host, length n_frames + 1), per observation its raw hkl as the sort key until the sort has used - // it, and its delta_phi as integrated until the smoothing has + // the resident ingest (IngestUpload .. IngestFinish): each frame's first reflection (host, length + // n_frames + 1); per reflection what the sort and the cut read - its raw hkl as the sort key, its + // frame position, resolution and background - until they have, and the place the cut gives it (as + // an exclusive prefix over the kept ones, n_keep in all) until the observations are taken from the + // reflections; per observation its delta_phi as integrated until the smoothing has used it std::vector in_src_start; CudaDevicePtr in_key; + CudaDevicePtr in_img, in_d, in_bkg; + CudaDevicePtr in_new_idx; + int in_n_keep = 0; CudaDevicePtr in_dphi; + // The reflections of the ingest's frames, through two pinned buffers of whole frames, each filled on + // all threads while the other is on its way to the card; take(cnt, base, src) takes cnt of them, + // the first being reflection base, apart on the card. + template + void StageReflections(const Reflection *const *refl, const int32_t *count, size_t nthreads, Take take); RotationScaleMergeGPU::IngestLayout in_layout; }; @@ -1362,7 +1378,6 @@ void RotationScaleMergeGPU::Impl::AllocObservations(int n_in) { zeta = Alloc(n); corr = Alloc(n); bkg = Alloc(n); var_bkg = Alloc(n); image_number = Alloc(n); d_obs = Alloc(n); - px_obs = Alloc(n); py_obs = Alloc(n); frame = Alloc(n); on_ice = Alloc(n); clipped = Alloc(n); @@ -1386,10 +1401,10 @@ void RotationScaleMergeGPU::Impl::AllocFrames(const int32_t *fstart, const int32 void RotationScaleMergeGPU::Impl::ReleaseFulls() { n_fulls = 0; - f_h = {}; f_k = {}; f_l = {}; f_frame = {}; f_group = {}; - f_I = {}; f_sigma = {}; f_d = {}; f_img = {}; f_px = {}; f_py = {}; + f_run = {}; f_peak = {}; f_frame = {}; f_group = {}; + f_I = {}; f_sigma = {}; f_d = {}; f_img = {}; f_var_bkg = {}; f_var_per_I = {}; f_capture = {}; f_on_ice = {}; f_clipped = {}; - f_corr = {}; f_partiality = {}; f_zeta = {}; f_obs = {}; + f_corr = {}; f_zeta = {}; f_obs = {}; f_frame_perm = {}; f_frame_start = {}; f_frame_count = {}; f_gperm = {}; f_gstart = {}; f_gcount = {}; m_em_mean = {}; m_hand = {}; m_has_hands = {}; m_rejected = {}; @@ -1434,8 +1449,6 @@ void RotationScaleMergeGPU::SetObsField(ObsField f, int offset, int count, const case ObsField::VarBkg: UploadChunk(d.var_bkg, offset, count, v, impl_->s()); break; case ObsField::ImageNumber: UploadChunk(d.image_number, offset, count, v, impl_->s()); break; case ObsField::D: UploadChunk(d.d_obs, offset, count, v, impl_->s()); break; - case ObsField::Px: UploadChunk(d.px_obs, offset, count, v, impl_->s()); break; - case ObsField::Py: UploadChunk(d.py_obs, offset, count, v, impl_->s()); break; } } @@ -1457,7 +1470,7 @@ void RotationScaleMergeGPU::CopyObservationsFrom(const RotationScaleMergeGPU &fr const size_t n = d.n_obs; copy(d.obs, f.obs, n); copy(d.zeta, f.zeta, n); copy(d.corr, f.corr, n); copy(d.bkg, f.bkg, n); copy(d.var_bkg, f.var_bkg, n); copy(d.image_number, f.image_number, n); - copy(d.d_obs, f.d_obs, n); copy(d.px_obs, f.px_obs, n); copy(d.py_obs, f.py_obs, n); + copy(d.d_obs, f.d_obs, n); copy(d.frame, f.frame, n); copy(d.on_ice, f.on_ice, n); copy(d.clipped, f.clipped, n); const size_t nr = std::max(1, d.n_runs); d.perm = d.Alloc(std::max(1, d.n_perm)); @@ -1619,7 +1632,7 @@ void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partialit MergeParams p{}; p.n_groups = ng; p.min_partiality = min_partiality; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); @@ -1702,7 +1715,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ p.reject_down = -std::expm1(-reject_nsigma * reject_mult_sd); p.reject_median = median.get(); p.reject_var_add = var_add.get(); - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.half = d_half.get(); @@ -1760,7 +1773,7 @@ void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, d MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.merged_I = d_merged_I.get(); @@ -1797,7 +1810,7 @@ void RotationScaleMergeGPU::MergeBijvoetExcess(double error_model_a, double erro MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); @@ -1829,7 +1842,7 @@ void RotationScaleMergeGPU::MergeGroupScatter(double *sum, double *sum_sq, doubl MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); @@ -1856,15 +1869,15 @@ void RotationScaleMergeGPU::MergeHalves(uint64_t *obs_key, uint8_t *half) { MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); const int obs_blocks = std::min(65535, (nf + BLK - 1) / BLK); const int grp_blocks = std::min(65535, (ng + BLK - 1) / BLK); - ObservationKeyKernel<<s()>>>(nf, p, d.f_h.get(), d.f_k.get(), d.f_l.get(), key.get(), - usable.get()); + ObservationKeyKernel<<s()>>>(nf, p, d.f_run.get(), d.rr_h.get(), d.rr_k.get(), + d.rr_l.get(), key.get(), usable.get()); CudaCheck(cudaGetLastError(), "observation key launch"); CudaCheck(cudaMemsetAsync(h.get(), 0, size_t(nf), impl_->s()), "memset halves"); HalvesKernel<<s()>>>(p, key.get(), usable.get(), h.get()); @@ -1981,7 +1994,6 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.max_frame_gap = max_frame_gap; p.obs = d.obs.get(); p.corr = d.corr.get(); p.bkg = d.bkg.get(); p.var_bkg = d.var_bkg.get(); p.image_number = d.image_number.get(); p.d = d.d_obs.get(); - p.px = d.px_obs.get(); p.py = d.py_obs.get(); p.frame = d.frame.get(); p.on_ice = d.on_ice.get(); p.clipped = d.clipped.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 = rr_group.get(); @@ -2011,23 +2023,21 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti // Allocate the fulls SoA and emit. const int nf = std::max(1, d.n_fulls); - d.f_h = d.Alloc(nf); d.f_k = d.Alloc(nf); d.f_l = d.Alloc(nf); + d.f_run = d.Alloc(nf); d.f_peak = d.Alloc(nf); d.f_frame = d.Alloc(nf); d.f_group = d.Alloc(nf); d.f_I = d.Alloc(nf); d.f_sigma = d.Alloc(nf); d.f_d = d.Alloc(nf); d.f_img = d.Alloc(nf); - d.f_px = d.Alloc(nf); d.f_py = d.Alloc(nf); d.f_var_bkg = d.Alloc(nf); d.f_var_per_I = d.Alloc(nf); d.f_capture = d.Alloc(nf); d.f_on_ice = d.Alloc(nf); d.f_clipped = d.Alloc(nf); - d.f_corr = d.Alloc(nf); d.f_partiality = d.Alloc(nf); + d.f_corr = d.Alloc(nf); d.f_zeta = d.Alloc(nf); d.f_obs = d.Alloc(nf); d.Upload(rr_offset, offset.data(), d.n_runs); p.rr_offset = 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_run = d.f_run.get(); p.f_peak = d.f_peak.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_px = d.f_px.get(); p.f_py = d.f_py.get(); p.f_var_bkg = d.f_var_bkg.get(); p.f_var_per_I = d.f_var_per_I.get(); p.f_capture = d.f_capture.get(); p.f_on_ice = d.f_on_ice.get(); p.f_clipped = d.f_clipped.get(); if (d.n_fulls > 0) { @@ -2038,7 +2048,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti return d.n_fulls; } -void RotationScaleMergeGPU::GetFulls(int32_t *h, int32_t *k, int32_t *l, float *I, float *sigma, float *d, +void RotationScaleMergeGPU::GetFulls(int32_t *run, float *I, float *sigma, float *d, float *image_number, int32_t *frame, uint8_t *on_ice, uint8_t *clipped, int32_t *group) const { DeviceGuard guard(impl_->device, impl_->available); @@ -2048,8 +2058,7 @@ void RotationScaleMergeGPU::GetFulls(int32_t *h, int32_t *k, int32_t *l, float * auto dl = [&](void *dst, const void *src, size_t bytes) { CopyAndWait(dst, src, bytes, cudaMemcpyDeviceToHost, impl_->s(), "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(run, dd.f_run.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)); @@ -2066,8 +2075,6 @@ void RotationScaleMergeGPU::ResetFullsScale() { // Unity model: partiality/prescaling_corr/zeta = 1 so coeff = mean; corr starts at 1. FillKernel<<s()>>>(d.f_corr.get(), nf, 1.0f); CudaCheck(cudaGetLastError(), "FillKernel launch"); - FillKernel<<s()>>>(d.f_partiality.get(), nf, 1.0f); - CudaCheck(cudaGetLastError(), "FillKernel launch"); FullsRecordKernel<<s()>>>(nf, d.f_I.get(), d.f_sigma.get(), d.f_obs.get()); CudaCheck(cudaGetLastError(), "FullsRecordKernel launch"); FillKernel<<s()>>>(d.f_zeta.get(), nf, 1.0f); @@ -2105,13 +2112,12 @@ void RotationScaleMergeGPU::GetFullsCorr(float *corr) const { cudaMemcpyDeviceToHost, impl_->s(), "download f_corr"); } -void RotationScaleMergeGPU::GetFullsPxPy(float *px, float *py) const { +void RotationScaleMergeGPU::GetFullsPeak(int32_t *peak) const { DeviceGuard guard(impl_->device, impl_->available); const auto &d = *impl_; if (d.n_fulls == 0) return; - const size_t bytes = size_t(d.n_fulls) * sizeof(float); - CopyAndWait(px, d.f_px.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_px"); - CopyAndWait(py, d.f_py.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_py"); + CopyAndWait(peak, d.f_peak.get(), size_t(d.n_fulls) * sizeof(int32_t), cudaMemcpyDeviceToHost, impl_->s(), + "download f_peak"); } void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I, float *capture) const { @@ -2128,7 +2134,7 @@ void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I, f void RotationScaleMergeGPU::FullsDownloaded() { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; - d.f_img = {}; d.f_px = {}; d.f_py = {}; d.f_clipped = {}; + d.f_img = {}; d.f_peak = {}; d.f_clipped = {}; d.f_obs = {}; d.f_zeta = {}; d.f_frame_perm = {}; d.f_frame_start = {}; d.f_frame_count = {}; d.group_mean = {}; } @@ -2456,19 +2462,10 @@ namespace { return (u & 0x80000000u) ? ~u : (u | 0x80000000u); } - struct IngestArrays { - ObsRecord *obs; - float *zeta, *corr, *bkg, *var_bkg, *image_number, *d, *px, *py, *dphi; - int32_t *frame; - uint8_t *on_ice, *clipped; - uint64_t *key; - }; - // One block per BLK reflections: the block reads its reflections as whole words (coalesced) into - // shared memory and each thread takes its own apart. The deterministic correction is the product of - // its three factors, in the host's order; the frame's flux is divided out once it is known. - __global__ void IngestExtractKernel(int n, int base, const Reflection *__restrict__ refl, IngestArrays a) { - __shared__ uint32_t buf[BLK * REFL_WORDS]; + // shared memory, and each thread then has its own at hand. + __device__ __forceinline__ const Reflection *StageBlock(int n, const Reflection *__restrict__ refl, + uint32_t *buf) { const int first = blockIdx.x * BLK; const int cnt = min(BLK, n - first); const uint32_t *src = reinterpret_cast(refl + first); @@ -2476,9 +2473,48 @@ namespace { buf[w] = src[w]; __syncthreads(); if (int(threadIdx.x) >= cnt) + return nullptr; + return reinterpret_cast(buf + threadIdx.x * REFL_WORDS); + } + + // What the sort and the cut read of every reflection. + struct IngestKeys { + float *image_number, *d, *bkg; + uint64_t *key; + }; + __global__ void IngestKeyKernel(int n, int base, const Reflection *__restrict__ refl, IngestKeys a) { + __shared__ uint32_t buf[BLK * REFL_WORDS]; + const Reflection *r = StageBlock(n, refl, buf); + if (!r) + return; + const int i = base + blockIdx.x * BLK + threadIdx.x; + a.image_number[i] = r->image_number; + a.d[i] = r->d; + a.bkg[i] = r->bkg; + a.key[i] = RawHklKey(r->h, r->k, r->l); + } + + struct IngestArrays { + ObsRecord *obs; + float *zeta, *corr, *bkg, *var_bkg, *image_number, *d, *dphi; + uint8_t *on_ice, *clipped; + const int32_t *new_idx; // the cut's place of each reflection (an exclusive prefix over the kept) + int n, n_keep; // reflections and observations + }; + + // The observations, each at the place the cut gave its reflection (a reflection the cut dropped has + // the place of the next kept one). The deterministic correction is the product of its three + // factors, in the host's order; the frame's flux is divided out once the frames are laid out. + __global__ void IngestExtractKernel(int cnt, int base, const Reflection *__restrict__ refl, IngestArrays a) { + __shared__ uint32_t buf[BLK * REFL_WORDS]; + const Reflection *rp = StageBlock(cnt, refl, buf); + if (!rp) + return; + const Reflection &r = *rp; + const int src = base + blockIdx.x * BLK + threadIdx.x; + const int i = a.new_idx[src]; + if ((src + 1 < a.n ? a.new_idx[src + 1] : a.n_keep) == i) return; - const Reflection &r = *reinterpret_cast(buf + threadIdx.x * REFL_WORDS); - const int i = base + first + threadIdx.x; a.obs[i] = ObsRecord{r.I, r.sigma, r.partiality, __fmul_rn(__fmul_rn(r.prescaling_corr, r.qe_corr), r.flight_corr)}; a.zeta[i] = r.zeta; @@ -2487,12 +2523,9 @@ namespace { a.var_bkg[i] = r.var_bkg; a.image_number[i] = r.image_number; a.d[i] = r.d; - a.px[i] = r.predicted_x; - a.py[i] = r.predicted_y; a.dphi[i] = r.delta_phi_deg; a.on_ice[i] = r.on_ice_ring ? 1 : 0; a.clipped[i] = (r.clipped ? 1 : 0) | (r.overloaded ? 2 : 0); // OBS_CLIPPED | OBS_OVERLOADED - a.key[i] = RawHklKey(r.h, r.k, r.l); } __global__ void IngestFrameKernel(int n_frames, const int32_t *__restrict__ start, int32_t *__restrict__ frame) { @@ -2606,12 +2639,6 @@ namespace { kept[r] = keep_run[r] ? count[r] : 0; } - template - __global__ void IngestCompactKernel(int n, const int32_t *__restrict__ keep, const int32_t *__restrict__ new_idx, - const T *__restrict__ in, T *__restrict__ out) { - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += gridDim.x * blockDim.x) - if (keep[i]) out[new_idx[i]] = in[i]; - } int Blocks(int n) { return std::max(1, std::min(65535, (n + BLK - 1) / BLK)); } @@ -2796,83 +2823,85 @@ namespace { } } -void RotationScaleMergeGPU::IngestUpload(int n_frames, const Reflection *const *refl, const int32_t *count, - size_t nthreads) { - DeviceGuard guard(impl_->device, impl_->available); - auto &d = *impl_; - d.n_frames = n_frames; - d.in_src_start.assign(n_frames + 1, 0); +template +void RotationScaleMergeGPU::Impl::StageReflections(const Reflection *const *refl, const int32_t *count, + size_t nthreads, Take take) { int max_count = 0; - for (int o = 0; o < n_frames; ++o) { - d.in_src_start[o + 1] = d.in_src_start[o] + count[o]; + for (int o = 0; o < n_frames; ++o) max_count = std::max(max_count, count[o]); - } - const int n = d.in_src_start[n_frames]; - d.n_obs = n; - size_t free_bytes = 0, total_bytes = 0; - CudaCheck(cudaMemGetInfo(&free_bytes, &total_bytes), "memory info"); - if (size_t(n) * OBS_BYTES > total_bytes) - ThrowTooLargeForGPU(n); - d.AllocObservations(n); - d.in_key = d.Alloc(std::max(1, n)); - d.in_dphi = d.Alloc(std::max(1, n)); - - CudaDevicePtr src_start; - d.Upload(src_start, d.in_src_start.data(), n_frames + 1); - IngestFrameKernel<<>>(n_frames, src_start.get(), - d.frame.get()); - CudaCheck(cudaGetLastError(), "IngestFrameKernel launch"); - - // Through two pinned buffers of whole frames, each filled on all threads while the other is on - // its way to the card and taken apart there. Where the card reads host memory directly - // (gpu_direct_host_read()), the kernel takes the pinned buffer apart in place, without the upload. + // Where the card reads host memory directly (gpu_direct_host_read()), the kernel takes the pinned + // buffer apart in place, without the upload. constexpr size_t STAGE_BYTES = 64 << 20; const bool direct = gpu_direct_host_read(); const int cap = std::max(static_cast(STAGE_BYTES / sizeof(Reflection)), max_count); CudaHostPtr stage[2] = {CudaHostPtr(cap), CudaHostPtr(cap)}; CudaDevicePtr staged[2]; if (!direct) { - staged[0] = d.Alloc(cap); - staged[1] = d.Alloc(cap); + staged[0] = Alloc(cap); + staged[1] = Alloc(cap); } CudaEvent copied[2]; bool in_flight[2] = {false, false}; - IngestArrays a{d.obs.get(), d.zeta.get(), d.corr.get(), d.bkg.get(), d.var_bkg.get(), d.image_number.get(), - d.d_obs.get(), d.px_obs.get(), d.py_obs.get(), d.in_dphi.get(), d.frame.get(), - d.on_ice.get(), d.clipped.get(), d.in_key.get()}; int b = 0; for (int f0 = 0; f0 < n_frames; ) { int f1 = f0 + 1; - while (f1 < n_frames && d.in_src_start[f1 + 1] - d.in_src_start[f0] <= cap) + while (f1 < n_frames && in_src_start[f1 + 1] - in_src_start[f0] <= cap) ++f1; - const int base = d.in_src_start[f0]; - const int cnt = d.in_src_start[f1] - base; + const int base = in_src_start[f0]; + const int cnt = in_src_start[f1] - base; if (in_flight[b]) CudaCheck(cudaEventSynchronize(copied[b]), "staging buffer"); Reflection *dst = stage[b].get(); ParallelFor(f1 - f0, nthreads, [&](int fo) { const int o = f0 + fo; - std::copy(refl[o], refl[o] + count[o], dst + (d.in_src_start[o] - base)); + std::copy(refl[o], refl[o] + count[o], dst + (in_src_start[o] - base)); }); if (cnt > 0) { const Reflection *src = dst; if (!direct) { CudaCheck(cudaMemcpyAsync(staged[b].get(), dst, size_t(cnt) * sizeof(Reflection), - cudaMemcpyHostToDevice, d.s()), "upload reflections"); - CudaCheck(cudaEventRecord(copied[b], d.s()), "upload reflections"); + cudaMemcpyHostToDevice, s()), "upload reflections"); + CudaCheck(cudaEventRecord(copied[b], s()), "upload reflections"); src = staged[b].get(); } in_flight[b] = true; - IngestExtractKernel<<<(cnt + BLK - 1) / BLK, BLK, 0, d.s()>>>(cnt, base, src, a); - CudaCheck(cudaGetLastError(), "IngestExtractKernel launch"); + take(cnt, base, src); // Read in place, the buffer is free again only once the kernel has run. if (direct) - CudaCheck(cudaEventRecord(copied[b], d.s()), "extract reflections"); + CudaCheck(cudaEventRecord(copied[b], s()), "extract reflections"); } b ^= 1; f0 = f1; } - CudaCheck(cudaStreamSynchronize(d.s()), "upload reflections"); + CudaCheck(cudaStreamSynchronize(s()), "upload reflections"); +} + +void RotationScaleMergeGPU::IngestUpload(int n_frames, const Reflection *const *refl, const int32_t *count, + size_t nthreads) { + DeviceGuard guard(impl_->device, impl_->available); + auto &d = *impl_; + d.n_frames = n_frames; + d.in_src_start.assign(n_frames + 1, 0); + for (int o = 0; o < n_frames; ++o) + d.in_src_start[o + 1] = d.in_src_start[o] + count[o]; + const int n = d.in_src_start[n_frames]; + d.n_obs = n; + size_t free_bytes = 0, total_bytes = 0; + CudaCheck(cudaMemGetInfo(&free_bytes, &total_bytes), "memory info"); + if (size_t(n) * OBS_BYTES > total_bytes) + ThrowTooLargeForGPU(n); + // See AllocObservations: these are the first of the scaling's buffers. + cuda_trim_allocation_pool(); + // Only what the sort and the cut read; the observations are taken from the reflections once the cut + // has decided which they are (IngestObservations), so the two are never on the card together. + const int na = std::max(1, n); + d.in_img = d.Alloc(na); d.in_d = d.Alloc(na); d.in_bkg = d.Alloc(na); + d.in_key = d.Alloc(na); + const IngestKeys a{d.in_img.get(), d.in_d.get(), d.in_bkg.get(), d.in_key.get()}; + d.StageReflections(refl, count, nthreads, [&](int cnt, int base, const Reflection *src) { + IngestKeyKernel<<<(cnt + BLK - 1) / BLK, BLK, 0, d.s()>>>(cnt, base, src, a); + CudaCheck(cudaGetLastError(), "IngestKeyKernel launch"); + }); } void RotationScaleMergeGPU::IngestMeanBackground(double *mean_bkg) { @@ -2881,26 +2910,19 @@ void RotationScaleMergeGPU::IngestMeanBackground(double *mean_bkg) { CudaDevicePtr start; d.Upload(start, d.in_src_start.data(), d.n_frames + 1); CudaDevicePtr mean = d.Alloc(std::max(1, d.n_frames)); - IngestMeanBkgKernel<<>>(d.n_frames, start.get(), d.bkg.get(), mean.get()); + IngestMeanBkgKernel<<>>(d.n_frames, start.get(), d.in_bkg.get(), mean.get()); CudaCheck(cudaGetLastError(), "IngestMeanBkgKernel launch"); if (d.n_frames > 0) CopyAndWait(mean_bkg, mean.get(), size_t(d.n_frames) * sizeof(double), cudaMemcpyDeviceToHost, d.s(), "download mean background"); + d.in_bkg = {}; } -int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional d_min, - std::optional d_max) { +int RotationScaleMergeGPU::IngestSort(std::optional d_min, std::optional d_max) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; const int n = d.n_obs; const int na = std::max(1, n); - { - CudaDevicePtr flux; - d.Upload(flux, frame_flux, d.n_frames); - IngestFluxKernel<<>>(n, d.frame.get(), flux.get(), d.obs.get()); - CudaCheck(cudaGetLastError(), "IngestFluxKernel launch"); - CudaCheck(cudaStreamSynchronize(d.s()), "flux"); - } // (h, k, l, image_number, index) is a total order, so two stable radix passes - by the frame position // with the index as the value, then by the raw hkl - give the one permutation std::sort gives on the @@ -2926,10 +2948,11 @@ int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional ikey = d.Alloc(na), ikey_sorted = d.Alloc(na); CudaDevicePtr idx = d.Alloc(na); - IngestImageKeyKernel<<>>(n, d.image_number.get(), ikey.get(), idx.get()); + IngestImageKeyKernel<<>>(n, d.in_img.get(), ikey.get(), idx.get()); CudaCheck(cudaGetLastError(), "IngestImageKeyKernel launch"); if (n > 0) sort_pairs(ikey, ikey_sorted, idx, by_image, 32); } + d.in_img = {}; CudaDevicePtr key = d.Alloc(na); IngestGatherKeyKernel<<>>(n, by_image.get(), d.in_key.get(), key.get()); CudaCheck(cudaGetLastError(), "IngestGatherKeyKernel launch"); @@ -2984,7 +3007,7 @@ int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional keep_run = d.Alloc(nra); const bool limited = d_min.has_value() || d_max.has_value(); IngestRunKernel<<>>(n_runs, run_start.get(), run_count.get(), perm.get(), - d.d_obs.get(), d_min.has_value(), d_min.value_or(0.0), + d.in_d.get(), d_min.has_value(), d_min.value_or(0.0), d_max.has_value(), d_max.value_or(0.0), run_d.get(), keep_run.get()); CudaCheck(cudaGetLastError(), "IngestRunKernel launch"); @@ -3043,6 +3066,7 @@ int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional; - CudaDevicePtr out = d.Alloc(nka); - IngestCompactKernel<<>>(n, keep.get(), new_idx.get(), field.get(), out.get()); - CudaCheck(cudaGetLastError(), "IngestCompactKernel launch"); - CudaCheck(cudaStreamSynchronize(d.s()), "cut"); - field = std::move(out); - }; - cut(d.obs); cut(d.zeta); cut(d.corr); cut(d.bkg); cut(d.var_bkg); cut(d.image_number); cut(d.d_obs); - cut(d.px_obs); cut(d.py_obs); cut(d.frame); cut(d.on_ice); cut(d.clipped); cut(d.in_dphi); - } + // Where each observation goes, for IngestObservations to take them from the reflections. + d.in_new_idx = std::move(new_idx); + d.in_n_keep = n_keep; d.n_obs = n_keep; d.AllocFrames(L.frame_start.data(), L.frame_count.data()); @@ -3094,6 +3108,51 @@ void RotationScaleMergeGPU::IngestGetLayout(IngestLayout &out) const { out = impl_->in_layout; } +void RotationScaleMergeGPU::IngestObservations(const Reflection *const *refl, const int32_t *count, size_t nthreads, + const double *frame_flux, float *px, float *py) { + DeviceGuard guard(impl_->device, impl_->available); + auto &d = *impl_; + const int n_keep = d.in_n_keep; + d.AllocObservations(n_keep); + d.in_dphi = d.Alloc(std::max(1, n_keep)); + const IngestArrays a{d.obs.get(), d.zeta.get(), d.corr.get(), d.bkg.get(), d.var_bkg.get(), d.image_number.get(), + d.d_obs.get(), d.in_dphi.get(), d.on_ice.get(), + d.clipped.get(), d.in_new_idx.get(), d.in_src_start[d.n_frames], n_keep}; + d.StageReflections(refl, count, nthreads, [&](int cnt, int base, const Reflection *src) { + IngestExtractKernel<<<(cnt + BLK - 1) / BLK, BLK, 0, d.s()>>>(cnt, base, src, a); + CudaCheck(cudaGetLastError(), "IngestExtractKernel launch"); + }); + // The predicted positions, which only the fulls' and only on the host read, stay there. + const int n = d.in_src_start[d.n_frames]; + std::vector new_idx(n); + if (n > 0) + CopyAndWait(new_idx.data(), d.in_new_idx.get(), size_t(n) * sizeof(int32_t), cudaMemcpyDeviceToHost, d.s(), + "download places"); + d.in_new_idx = {}; + ParallelFor(d.n_frames, nthreads, [&](int o) { + for (int j = 0; j < count[o]; ++j) { + const int src = d.in_src_start[o] + j, i = new_idx[src]; + if ((src + 1 < n ? new_idx[src + 1] : n_keep) == i) continue; + px[i] = refl[o][j].predicted_x; + py[i] = refl[o][j].predicted_y; + } + }); + + // The kept observations are in frame order, as the reflections were. + std::vector start(d.n_frames + 1, n_keep); + for (int o = 0; o < d.n_frames; ++o) start[o] = d.in_layout.frame_start[o]; + CudaDevicePtr d_start; + d.Upload(d_start, start.data(), d.n_frames + 1); + IngestFrameKernel<<>>(d.n_frames, d_start.get(), + d.frame.get()); + CudaCheck(cudaGetLastError(), "IngestFrameKernel launch"); + CudaDevicePtr flux; + d.Upload(flux, frame_flux, d.n_frames); + IngestFluxKernel<<>>(n_keep, d.frame.get(), flux.get(), d.obs.get()); + CudaCheck(cudaGetLastError(), "IngestFluxKernel launch"); + CudaCheck(cudaStreamSynchronize(d.s()), "flux"); +} + int64_t RotationScaleMergeGPU::IngestDeltaPhi(const DeltaPhiFrame &g) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; @@ -3173,14 +3232,19 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector &key_in, CudaDevicePtr &key_out, + CudaDevicePtr &val_in, CudaDevicePtr &val_out) { + if (nf > 0) IotaKernel<<<(nf + BLK - 1) / BLK, BLK, 0, d.s()>>>(nf, val_in.get()); + CudaCheck(cudaGetLastError(), "IotaKernel launch"); + cub::DoubleBuffer keys(key_in.get(), key_out.get()), values(val_in.get(), val_out.get()); size_t bytes = 0; - CudaCheck(cub::DeviceRadixSort::SortPairs(nullptr, bytes, key_in, key_out, val_in, val_out, nf, 0, 32, d.s()), - "sort temp size"); + CudaCheck(cub::DeviceRadixSort::SortPairs(nullptr, bytes, keys, values, nf, 0, 32, d.s()), "sort temp size"); CudaDevicePtr temp = d.Alloc(std::max(1, bytes)); - CudaCheck(cub::DeviceRadixSort::SortPairs(temp.get(), bytes, key_in, key_out, val_in, val_out, nf, 0, 32, d.s()), - "sort"); + CudaCheck(cub::DeviceRadixSort::SortPairs(temp.get(), bytes, keys, values, nf, 0, 32, d.s()), "sort"); CudaCheck(cudaStreamSynchronize(d.s()), "sort"); + if (values.Current() != val_out.get()) std::swap(val_in, val_out); }; // The counts come back for the starts, a prefix the host forms in order. const auto count = [&](const int32_t *bin, int n_bins, CudaDevicePtr &out, std::vector &host) { @@ -3199,8 +3263,6 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector idx = d.Alloc(nfa), key = d.Alloc(nfa), sorted = d.Alloc(nfa); - if (nf > 0) IotaKernel<<<(nf + BLK - 1) / BLK, BLK, 0, d.s()>>>(nf, idx.get()); - CudaCheck(cudaGetLastError(), "IotaKernel launch"); std::vector fcount; count(d.f_frame.get(), nfr, d.f_frame_count, fcount); @@ -3208,7 +3270,11 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector fstart = prefix(fcount); d.Upload(d.f_frame_start, fstart.data(), nfr); d.f_frame_perm = d.Alloc(nfa); - if (nf > 0) stable_sort(d.f_frame.get(), key.get(), idx.get(), d.f_frame_perm.get()); + if (nf > 0) { + CopyAndWait(key.get(), d.f_frame.get(), size_t(nf) * sizeof(int32_t), cudaMemcpyDeviceToDevice, d.s(), + "frame key"); + stable_sort(key, sorted, idx, d.f_frame_perm); + } std::vector gcount; count(d.f_group.get(), ng, d.f_gcount, gcount); @@ -3218,7 +3284,7 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector>>(nf, d.f_group.get(), key.get()); CudaCheck(cudaGetLastError(), "FullGroupKeyKernel launch"); CudaDevicePtr by_group = d.Alloc(nfa); - if (nf > 0) stable_sort(key.get(), sorted.get(), idx.get(), by_group.get()); + if (nf > 0) stable_sort(key, sorted, idx, by_group); d.f_gperm = d.Alloc(std::max(1, n_grouped)); gperm.resize(n_grouped); if (n_grouped > 0) { diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index b5d0c64b8..054439b13 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -40,12 +40,12 @@ public: // One slice ([offset, offset+count) of n_obs) of one immutable per-observation field, so the // host can stage a bounded chunk at a time out of its AoS record: staged whole and all at once, // the fourteen arrays are more than half the entire observation payload held live again. - // Corr0 seeds the resident, mutable corr (SetCorr refreshes it). Bkg/VarBkg/ImageNumber/D/Px/Py + // Corr0 seeds the resident, mutable corr (SetCorr refreshes it). Bkg/VarBkg/ImageNumber/D // are the combine inputs (image-local background, fractional frame position for event - // contiguity, resolution, and the predicted detector position carried to the full for the - // absorption surface). + // contiguity and resolution; the predicted detector position carried to the full for the + // absorption surface stays on the host, GetFullsPeak). enum class ObsField { I, Sigma, PrescalingCorr, Partiality, Zeta, Corr0, - Bkg, VarBkg, ImageNumber, D, Px, Py }; + Bkg, VarBkg, ImageNumber, D }; void SetObsField(ObsField f, int offset, int count, const float *v); void SetObsFrame(int offset, int count, const int32_t *frame); void SetObsOnIce(int offset, int count, const uint8_t *on_ice); @@ -60,19 +60,20 @@ public: // --- the resident ingest (RotationScaleMerge::IngestDevice with no observation array on the host) --- // The steps below build, from the reflections, exactly the observations and raw-hkl runs that // RotationScaleMerge::BuildInRangeObservations builds on the host, and leave them where the uploads - // above and SetRawRuns would: the reflections cross the bus once, and are sorted and cut on the card. + // above and SetRawRuns would: the reflections are sorted and cut on the card, and cross the bus + // twice - once for what the sort and the cut read, once for the observations the cut kept. // The reflections of n_frames frames, frame o's count[o] of them at refl[o], uploaded in frame order; - // every field the ingest and the merge read is kept per reflection. + // what the sort, the cut and the background below read is kept per reflection. void IngestUpload(int n_frames, const Reflection *const *refl, const int32_t *count, size_t nthreads); // Each frame's mean background over its reflections with a finite, positive one, summed in their // order (NaN where there is none): the incident-flux meter. Length n_frames. void IngestMeanBackground(double *mean_bkg); - // Divide every reflection's deterministic correction by its frame's flux (length n_frames); order the - // reflections by (raw h,k,l, image_number, index) and cut that order into raw-hkl runs, each with the - // resolution of its first member that has a usable one; and drop the runs outside [d_min, d_max] - // (either may be absent) with their reflections. Returns how many observations are left. - int IngestSort(const double *frame_flux, std::optional d_min, std::optional d_max); + // Order the reflections by (raw h,k,l, image_number, index) and cut that order into raw-hkl runs, + // each with the resolution of its first member that has a usable one; and drop the runs outside + // [d_min, d_max] (either may be absent) with their reflections. Returns how many observations are + // left. Nothing per observation is on the card yet (IngestObservations). + int IngestSort(std::optional d_min, std::optional d_max); // The layout IngestSort left: the reflection and run counts before the resolution cut, the // observations' frame CSR, and the runs that survived it. struct IngestLayout { @@ -83,6 +84,12 @@ public: std::vector run_d; }; void IngestGetLayout(IngestLayout &out) const; + // The observations IngestSort kept, taken from the same reflections (as IngestUpload had them), in + // its order, each with its deterministic correction divided by its frame's flux (length n_frames). + // px/py (length IngestSort's count) receive each observation's predicted position, which is not + // kept on the card (GetFullsPeak). + void IngestObservations(const Reflection *const *refl, const int32_t *count, size_t nthreads, + const double *frame_flux, float *px, float *py); // The ingest's geometry smoothing, on the observations IngestSort left (RotationScaleMerge::IngestOnDevice): // each delta_phi re-solved against its frame's smoothed lattice where the frame has one - recip holds // a*, b*, c* per frame (9 floats), ok which frames, and the goniometer frame as SmoothedGeometry @@ -229,15 +236,16 @@ public: int Combine(const int32_t *rawrun_group, double min_partiality, double capture_uncertainty_coeff, double min_captured_fraction, float max_frame_gap, int64_t &n_overloaded); - // Download the combined fulls SoA (length = Combine()'s return). The working corr is downloaded - // separately by GetFullsCorr (it is only meaningful after the fulls scaling; otherwise the caller sets it). - void GetFulls(int32_t *h, int32_t *k, int32_t *l, float *I, float *sigma, float *d, + // Download the combined fulls SoA (length = Combine()'s return); run is the raw-hkl run each full is + // of. The working corr is downloaded separately by GetFullsCorr (it is only meaningful after the + // fulls scaling; otherwise the caller sets it). + void GetFulls(int32_t *run, float *I, float *sigma, float *d, float *image_number, int32_t *frame, uint8_t *on_ice, uint8_t *clipped, int32_t *group) const; - // Download the fulls' predicted detector position (peak partial's px/py), for the host absorption - // surface. Length = n_fulls. - void GetFullsPxPy(float *px, float *py) const; + // Download the observation each full's predicted detector position is that of (its peak partial), + // for the host absorption surface; the positions themselves stay on the host. Length = n_fulls. + void GetFullsPeak(int32_t *peak) const; // Download the fulls' variance model, var(I) = var_bkg + var_per_I * I + (capture * I)^2. Length = n_fulls. void GetFullsVariance(float *var_bkg, float *var_per_I, float *capture) const; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index e39a262e7..e9ce40980 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1162,7 +1162,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // free: the answer must not depend on the moment. What else stays on the card meanwhile - the // context, the run's indexers - was measured at up to 2.3 GB; a speculative pass running beside // holds more, but a merge that does not fit beside it waits for it to end. On a 16 GB card that - // is the data sets up to about 70 M observations at high multiplicity; the large ones merge + // is the data sets up to about 90 M observations at high multiplicity; the large ones merge // their arms one after the other. if (!geometry_prepass && !experiment_.GetGemmiSpaceGroup().has_value() && rsm->GetSearchMinZeta() > 0.0 && cards > 0 && rsm->OnGPU()) { diff --git a/tests/RotationIngestGPUTest.cpp b/tests/RotationIngestGPUTest.cpp index bb78e8e48..7b8588480 100644 --- a/tests/RotationIngestGPUTest.cpp +++ b/tests/RotationIngestGPUTest.cpp @@ -97,7 +97,7 @@ TEST_CASE("RotationIngestGPU_SortAndCutMatchHost", "[RotationScale][gpu]") { std::vector flux(n_frames); for (int f = 0; f < n_frames; ++f) flux[f] = 0.9 + 0.01 * f; - const int n_keep = gpu.IngestSort(flux.data(), d_min, d_max); + const int n_keep = gpu.IngestSort(d_min, d_max); // The host ingest: one total order, its raw-hkl runs, each with its first usable resolution, and // the runs the range keeps with their observations. @@ -160,6 +160,14 @@ TEST_CASE("RotationIngestGPU_SortAndCutMatchHost", "[RotationScale][gpu]") { s += c; } + std::vector px(n_keep), py(n_keep); + gpu.IngestObservations(ptr.data(), count.data(), 4, flux.data(), px.data(), py.data()); + for (int i = 0; i < n; ++i) + if (keep[i]) { + CHECK(px[new_idx[i]] == all[i]->predicted_x); + CHECK(py[new_idx[i]] == all[i]->predicted_y); + } + // The rocking events' spans, walked over each run's members in the sorted order, and corr as // integrated in the observations' order (no frame has a mosaicity to recompute partiality from). constexpr float GAP = 1.0f;