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

Merge scale-footprint round 2 into rc175: two-pass scaling ingest with the GPU/CPU decision before the card is touched, px/py and h,k,l kept on the host, estimate re-fitted to 56 B per observation + 99 B per distinct hkl, output unchanged

Browse Source 
8qaw's scaling engine peak 10.76 -> 8.16 GB; p.hkl byte-identical on 15 battery sets.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-15 11:14:49 +02:00
2 parents 9d8b06525f f46dc303a3
commit bd7464d705
7 files changed
+346 -222
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/scale_merge/ErrorModelGPU.cu
+12 -9
View File
@@ -89,8 +89,8 @@ ErrorModelFit FitErrorModelGPU(const std::vector<ErrorModelSample> &pool) {
return fit;
const int n = static_cast<int>(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<ErrorModelSample> &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<int32_t> idx(n, ALLOC), idx_next(n, ALLOC);
CudaDevicePtr<double> key(n, ALLOC), key_sorted(n, ALLOC);
IotaKernel<<<Blocks(n), THREADS, 0, stream>>>(n, idx);
cuda_err(cudaGetLastError());
cub::DoubleBuffer<double> keys(key.get(), key_sorted.get());
cub::DoubleBuffer<int32_t> 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<uint8_t> scratch(std::max<size_t>(bytes, 1), ALLOC);
for (int field = 0; field < 4; ++field) {
KeyKernel<<<Blocks(n), THREADS, 0, stream>>>(n, d_pool, idx, field, key);
KeyKernel<<<Blocks(n), THREADS, 0, stream>>>(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<uint8_t>(); idx_next = CudaDevicePtr<int32_t>(); key = CudaDevicePtr<double>();
// After the last pass key_sorted holds I2/s2 in rank order.
CudaDevicePtr<double> s2(n, ALLOC), I2(n, ALLOC), dev2(n, ALLOC);
image_analysis/scale_merge/RotationScaleMerge.cpp
+68 -37
View File
@@ -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<const Reflection *> refl(n_frames);
std::vector<int32_t> 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<int32_t>(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<double> 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<ObsPosition>();
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<size_t>(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<size_t>(resident_ingest ? rocking_events_at_ingest : n_partials_obs);
const size_t bytes = BYTES_PER_OBS * static_cast<size_t>(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<float[]> fbuf[NFLOAT];
for (auto &b : fbuf) b.reset(new float[slice]);
std::unique_ptr<int32_t[]> 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<ObsPosition>();
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<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());
@@ -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.
image_analysis/scale_merge/RotationScaleMerge.h
+14 -6
View File
@@ -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<float> px, py; };
std::shared_ptr<const ObsPosition> 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<int32_t> 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<int32_t> h, k, l, frame, group;
std::vector<float> I, sigma, d, image_number, corr, px, py, var_bkg, var_per_I, capture;
std::vector<int32_t> run, peak, frame, group;
std::vector<float> I, sigma, d, image_number, corr, var_bkg, var_per_I, capture;
std::vector<uint8_t> 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
image_analysis/scale_merge/RotationScaleMergeGPU.cu
+217 -151
View File
@@ -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<float> bkg, var_bkg, image_number, d_obs, px_obs, py_obs;
CudaDevicePtr<float> bkg, var_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;
// 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, 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<int32_t> f_run, f_peak, f_frame, f_group;
CudaDevicePtr<float> f_I, f_sigma, f_d, f_img, f_var_bkg, f_var_per_I, f_capture;
CudaDevicePtr<uint8_t> 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<float> f_corr, f_partiality, f_zeta;
CudaDevicePtr<float> f_corr, f_zeta;
CudaDevicePtr<ObsRecord> f_obs;
CudaDevicePtr<int32_t> f_frame_perm, f_frame_start, f_frame_count;
CudaDevicePtr<int32_t> f_gperm, f_gstart, f_gcount;
@@ -1285,12 +1291,22 @@ struct RotationScaleMergeGPU::Impl {
CudaDevicePtr<int32_t> s_perm[3], s_seg_start[3];
CudaDevicePtr<double> s_A, s_sw, s_swI, s_tcross, s_tref2, s_cross, s_ref2;
CudaDevicePtr<double> 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<int32_t> in_src_start;
CudaDevicePtr<uint64_t> in_key;
CudaDevicePtr<float> in_img, in_d, in_bkg;
CudaDevicePtr<int32_t> in_new_idx;
int in_n_keep = 0;
CudaDevicePtr<float> 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 <typename Take>
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<float>(n); corr = Alloc<float>(n);
bkg = Alloc<float>(n); var_bkg = Alloc<float>(n);
image_number = Alloc<float>(n); d_obs = Alloc<float>(n);
px_obs = Alloc<float>(n); py_obs = Alloc<float>(n);
frame = Alloc<int32_t>(n);
on_ice = Alloc<uint8_t>(n);
clipped = Alloc<uint8_t>(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<int32_t>(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<<<obs_blocks, BLK, 0, impl_->s()>>>(nf, p, d.f_h.get(), d.f_k.get(), d.f_l.get(), key.get(),
usable.get());
ObservationKeyKernel<<<obs_blocks, BLK, 0, impl_->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<<<grp_blocks, BLK, 0, impl_->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<int32_t>(nf); d.f_k = d.Alloc<int32_t>(nf); d.f_l = d.Alloc<int32_t>(nf);
d.f_run = d.Alloc<int32_t>(nf); d.f_peak = d.Alloc<int32_t>(nf);
d.f_frame = d.Alloc<int32_t>(nf); d.f_group = d.Alloc<int32_t>(nf);
d.f_I = d.Alloc<float>(nf); d.f_sigma = d.Alloc<float>(nf);
d.f_d = d.Alloc<float>(nf); d.f_img = d.Alloc<float>(nf);
d.f_px = d.Alloc<float>(nf); d.f_py = d.Alloc<float>(nf);
d.f_var_bkg = d.Alloc<float>(nf); d.f_var_per_I = d.Alloc<float>(nf); d.f_capture = d.Alloc<float>(nf);
d.f_on_ice = d.Alloc<uint8_t>(nf); d.f_clipped = d.Alloc<uint8_t>(nf);
d.f_corr = d.Alloc<float>(nf); d.f_partiality = d.Alloc<float>(nf);
d.f_corr = d.Alloc<float>(nf);
d.f_zeta = d.Alloc<float>(nf);
d.f_obs = d.Alloc<ObsRecord>(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<<<obs_blocks, BLK, 0, impl_->s()>>>(d.f_corr.get(), nf, 1.0f);
CudaCheck(cudaGetLastError(), "FillKernel launch");
FillKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(d.f_partiality.get(), nf, 1.0f);
CudaCheck(cudaGetLastError(), "FillKernel launch");
FullsRecordKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(nf, d.f_I.get(), d.f_sigma.get(), d.f_obs.get());
CudaCheck(cudaGetLastError(), "FullsRecordKernel launch");
FillKernel<<<obs_blocks, BLK, 0, impl_->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<const uint32_t *>(refl + first);
@@ -2476,9 +2473,48 @@ namespace {
buf[w] = src[w];
__syncthreads();
if (int(threadIdx.x) >= cnt)
return nullptr;
return reinterpret_cast<const Reflection *>(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<const Reflection *>(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 <typename T>
__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 <typename Take>
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<uint64_t>(std::max(1, n));
d.in_dphi = d.Alloc<float>(std::max(1, n));
CudaDevicePtr<int32_t> src_start;
d.Upload(src_start, d.in_src_start.data(), n_frames + 1);
IngestFrameKernel<<<std::max(1, std::min(n_frames, 65535)), BLK, 0, d.s()>>>(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<int>(STAGE_BYTES / sizeof(Reflection)), max_count);
CudaHostPtr<Reflection> stage[2] = {CudaHostPtr<Reflection>(cap), CudaHostPtr<Reflection>(cap)};
CudaDevicePtr<Reflection> staged[2];
if (!direct) {
staged[0] = d.Alloc<Reflection>(cap);
staged[1] = d.Alloc<Reflection>(cap);
staged[0] = Alloc<Reflection>(cap);
staged[1] = Alloc<Reflection>(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<float>(na); d.in_d = d.Alloc<float>(na); d.in_bkg = d.Alloc<float>(na);
d.in_key = d.Alloc<uint64_t>(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<int32_t> start;
d.Upload(start, d.in_src_start.data(), d.n_frames + 1);
CudaDevicePtr<double> mean = d.Alloc<double>(std::max(1, d.n_frames));
IngestMeanBkgKernel<<<Blocks(d.n_frames), BLK, 0, d.s()>>>(d.n_frames, start.get(), d.bkg.get(), mean.get());
IngestMeanBkgKernel<<<Blocks(d.n_frames), BLK, 0, d.s()>>>(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<double> d_min,
std::optional<double> d_max) {
int RotationScaleMergeGPU::IngestSort(std::optional<double> d_min, std::optional<double> 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<double> flux;
d.Upload(flux, frame_flux, d.n_frames);
IngestFluxKernel<<<Blocks(n), BLK, 0, d.s()>>>(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<do
{
CudaDevicePtr<uint32_t> ikey = d.Alloc<uint32_t>(na), ikey_sorted = d.Alloc<uint32_t>(na);
CudaDevicePtr<int32_t> idx = d.Alloc<int32_t>(na);
IngestImageKeyKernel<<<Blocks(n), BLK, 0, d.s()>>>(n, d.image_number.get(), ikey.get(), idx.get());
IngestImageKeyKernel<<<Blocks(n), BLK, 0, d.s()>>>(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<uint64_t> key = d.Alloc<uint64_t>(na);
IngestGatherKeyKernel<<<Blocks(n), BLK, 0, d.s()>>>(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<do
CudaDevicePtr<int32_t> keep_run = d.Alloc<int32_t>(nra);
const bool limited = d_min.has_value() || d_max.has_value();
IngestRunKernel<<<Blocks(n_runs), BLK, 0, d.s()>>>(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<do
// The runs are cut: what was kept of them is in rr_*, so the uncut ones and the old order go before
// the observations are cut below, one field held twice at a time.
perm = {}; run_key = {}; run_count = {}; run_start = {}; run_d = {}; keep_run = {}; run_idx = {}; run_at = {};
d.in_d = {};
// The frame CSR of what survived.
auto &L = d.in_layout;
@@ -3063,19 +3087,9 @@ int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional<do
L.frame_start.assign(d.n_frames, 0);
for (int o = 0, at = 0; o < d.n_frames; ++o) { L.frame_start[o] = at; at += L.frame_count[o]; }
// The observations themselves, cut field by field, so that only one field is ever held twice.
if (n_keep < n) {
const auto cut = [&](auto &field) {
using T = std::remove_reference_t<decltype(*field.get())>;
CudaDevicePtr<T> out = d.Alloc<T>(nka);
IngestCompactKernel<T><<<Blocks(n), BLK, 0, d.s()>>>(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<float>(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<int32_t> 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<int32_t> 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<int32_t> d_start;
d.Upload(d_start, start.data(), d.n_frames + 1);
IngestFrameKernel<<<std::max(1, std::min(d.n_frames, 65535)), BLK, 0, d.s()>>>(d.n_frames, d_start.get(),
d.frame.get());
CudaCheck(cudaGetLastError(), "IngestFrameKernel launch");
CudaDevicePtr<double> flux;
d.Upload(flux, frame_flux, d.n_frames);
IngestFluxKernel<<<Blocks(n_keep), BLK, 0, d.s()>>>(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<int
auto &d = *impl_;
const int nf = d.n_fulls, nfr = d.n_frames, ng = d.n_groups;
const int nfa = std::max(1, nf);
const auto stable_sort = [&](const int32_t *key_in, int32_t *key_out, const int32_t *val_in, int32_t *val_out) {
// The index order of the keys, stably sorted, left in val_out: back and forth between the two
// buffers of each, as no input is needed after it, so that cub holds no copy of its own.
const auto stable_sort = [&](CudaDevicePtr<int32_t> &key_in, CudaDevicePtr<int32_t> &key_out,
CudaDevicePtr<int32_t> &val_in, CudaDevicePtr<int32_t> &val_out) {
if (nf > 0) IotaKernel<<<(nf + BLK - 1) / BLK, BLK, 0, d.s()>>>(nf, val_in.get());
CudaCheck(cudaGetLastError(), "IotaKernel launch");
cub::DoubleBuffer<int32_t> 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<uint8_t> temp = d.Alloc<uint8_t>(std::max<size_t>(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<int32_t> &out, std::vector<int32_t> &host) {
@@ -3199,8 +3263,6 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector<int
return start;
};
CudaDevicePtr<int32_t> idx = d.Alloc<int32_t>(nfa), key = d.Alloc<int32_t>(nfa), sorted = d.Alloc<int32_t>(nfa);
if (nf > 0) IotaKernel<<<(nf + BLK - 1) / BLK, BLK, 0, d.s()>>>(nf, idx.get());
CudaCheck(cudaGetLastError(), "IotaKernel launch");
std::vector<int32_t> 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<int
const std::vector<int32_t> fstart = prefix(fcount);
d.Upload(d.f_frame_start, fstart.data(), nfr);
d.f_frame_perm = d.Alloc<int32_t>(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<int32_t> gcount;
count(d.f_group.get(), ng, d.f_gcount, gcount);
@@ -3218,7 +3284,7 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector<int
FullGroupKeyKernel<<<Blocks(nf), BLK, 0, d.s()>>>(nf, d.f_group.get(), key.get());
CudaCheck(cudaGetLastError(), "FullGroupKeyKernel launch");
CudaDevicePtr<int32_t> by_group = d.Alloc<int32_t>(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<int32_t>(std::max(1, n_grouped));
gperm.resize(n_grouped);
if (n_grouped > 0) {
image_analysis/scale_merge/RotationScaleMergeGPU.h
+25 -17
View File
@@ -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<double> d_min, std::optional<double> 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<double> d_min, std::optional<double> 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<float> 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;
rugnux/RugnuxScaleMerge.cpp
+1 -1
View File
@@ -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()) {
tests/RotationIngestGPUTest.cpp
+9 -1
View File
@@ -97,7 +97,7 @@ TEST_CASE("RotationIngestGPU_SortAndCutMatchHost", "[RotationScale][gpu]") {
std::vector<double> 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<float> 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;