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

rugnux: make the P1 merges of the tail beside the critical path

Browse Source 
The tail of the canonical pass made five merges one after another. Two of
them read nothing the space-group search or the in-symmetry merge decides:
the all-observation arm of the search and the P1 cross-check. They are now
made on a second RotationScaleMerge engine, ingested beside the run's own
before any merge writes per-frame values back, and taken where they were
made before; where the run re-ingests (a reindex, a cell change) they are
made on the run's engine as before.

- RotationScaleMerge::SetWriteBackPerFrameScale(false) keeps a merge that
  is not the run's answer from writing G/CC/mosaicity onto the outcomes.
  The P1 cross-check no longer overwrites them, so _plot.txt's scale_G,
  cc_to_merge and cc_n now describe the merge that was written (in the
  determined group) instead of the P1 cross-check; the cross-check also
  no longer leaks its scaling iteration count into the report.
- RotationScaleMergeGPU runs on its own non-blocking stream instead of
  the legacy NULL stream, so the two engines (and a probe pass beside
  them) do not serialise at every launch and synchronisation.
- The anisotropy analysis (mostly ScaledObservations) runs beside tNCS,
  twinning and the other report-only analyses.

p.mtz, p_P1.mtz, p.cif, p.hkl and p_unmerged.mtz are byte-identical on
myob/cytc/thau (GPU and CPU builds). Tail on cytc GPU 8.1 -> ~6.8-7.7 s
under a loaded box.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB
This commit is contained in:
leonarski_f
2026-10-03 10:06:22 +02:00
co-authored by Claude Opus 5.5
parent 2c97e654ed
commit 871347b7ad
4 changed files with 271 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+2 -1
View File
@@ -5699,7 +5699,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
ReducePartialGroupMeans(n_groups, partial_mean);
ComputePerFrameCC(partial_mean, cc, cc_n);
}
FinalizePerFrameScale(cc, cc_n, partial_scaled);
if (write_back_per_frame_scale)
FinalizePerFrameScale(cc, cc_n, partial_scaled);
// The filters below remove observations by zeroing corr, which is what takes an observation out of
// the 3D combine, the merge and the error model alike (excluding them from the ASU grouping is NOT
image_analysis/scale_merge/RotationScaleMerge.h
+7
View File
@@ -119,6 +119,12 @@ public:
// compares nothing, asks for it to be left out.
Result Run(bool for_search, bool full_stats, bool measure_cc_before_corrections);
// Whether Run() writes the per-frame G / CC / mosaicity back onto the outcomes (on by default). Off
// for a merge that is not the run's answer - the P1 cross-check - so the per-image table and the
// unmerged MTZ describe the merge that was written, and so an engine merging beside another one
// does not write the outcomes they share.
void SetWriteBackPerFrameScale(bool on) { write_back_per_frame_scale = on; }
// Override the high-resolution cut for the next Run() - used to gate the de-novo P1 search pass at
// <I/sigma> >= 1 without cutting the final in-symmetry merge. Reset to the manual limit afterwards.
void SetDMinLimit(std::optional<double> d_min_A) { d_min_limit = d_min_A; }
@@ -413,6 +419,7 @@ private:
std::unique_ptr<RotationScaleMergeGPU> gpu_;
bool gpu_active_ = false;
#endif
bool write_back_per_frame_scale = true; // see SetWriteBackPerFrameScale
// --- helpers (each a flat pass; see the .cpp) ---
// Turn the per-frame mean background under the reflections (accumulated by the ingest fill loop) into
image_analysis/scale_merge/RotationScaleMergeGPU.cu
+132 -115
View File
@@ -18,13 +18,16 @@ namespace {
constexpr int BLK = 256;
constexpr int MIN_REFLECTIONS = 20;
// Every kernel and copy here is queued on the legacy NULL stream, so - as in BeamCenterFFTGPU -
// no buffer comes from the pool. A pooled buffer is freed with cudaFreeAsync on the thread's
// non-blocking allocation stream, which is not ordered after the NULL stream. Each entry point
// below waits for its own work before it returns, so no free here has yet overtaken a read, but
// that holds only by that convention, and not at all for a free while unwinding from a failed
// call; nor can compute-sanitizer --track-stream-ordered-races see it, and it reported every
// reassigned merge buffer as a use-after-free. cudaFree synchronises the device first.
// Every kernel and copy here is queued on the instance's own stream (Impl::stream), not on the
// legacy NULL stream: two merges made side by side - the run's and the one it makes ahead of time -
// and the image analysis of a probe pass beside them would otherwise wait for each other's work at
// every launch and every synchronisation. As in BeamCenterFFTGPU no buffer comes from the pool. A
// pooled buffer is freed with cudaFreeAsync on the thread's allocation stream, which is not
// ordered after this one. Each entry point below waits for its own work before it returns, so no
// free here has yet overtaken a read, but that holds only by that convention, and not at all for a
// free while unwinding from a failed call; nor can compute-sanitizer --track-stream-ordered-races
// see it, and it reported every reassigned merge buffer as a use-after-free. cudaFree
// synchronises the device first.
constexpr CudaAlloc ALLOC = CudaAlloc::Synchronous;
__device__ __forceinline__ double SafeInvD(double x, double fallback) {
@@ -636,6 +639,16 @@ namespace {
}
}
void CudaCheck(cudaError_t e, const char *what);
// A copy on the instance's stream, waited for - what cudaMemcpy on the NULL stream was, without also
// waiting for every other stream on the card.
void CopyAndWait(void *dst, const void *src, size_t bytes, cudaMemcpyKind kind, cudaStream_t s,
const char *what) {
CudaCheck(cudaMemcpyAsync(dst, src, bytes, kind, s), what);
CudaCheck(cudaStreamSynchronize(s), what);
}
void CudaCheck(cudaError_t e, const char *what) {
if (e != cudaSuccess)
throw JFJochException(JFJochExceptionCategory::GPUCUDAError,
@@ -703,6 +716,9 @@ namespace {
}
struct RotationScaleMergeGPU::Impl {
// First, so it goes last: the buffers below are freed before the stream their work ran on.
std::unique_ptr<CudaStream> stream;
cudaStream_t s() const { return stream->get(); }
int device = 0; // the GPU this instance's buffers live on
bool available = false;
int n_obs = 0, n_frames = 0, n_groups = 0;
@@ -736,7 +752,7 @@ struct RotationScaleMergeGPU::Impl {
void Upload(CudaDevicePtr<T> &dst, const T *src, int n) const {
dst = Alloc<T>(std::max(1, n));
if (n > 0)
CudaCheck(cudaMemcpy(dst.get(), src, size_t(n) * sizeof(T), cudaMemcpyHostToDevice), "upload");
CopyAndWait(dst.get(), src, size_t(n) * sizeof(T), cudaMemcpyHostToDevice, s(), "upload");
}
// immutable per-obs
@@ -833,6 +849,7 @@ RotationScaleMergeGPU::RotationScaleMergeGPU() : impl_(std::make_unique<Impl>())
// can put the caller's device back instead of leaving the thread moved.
impl_->device = 0;
DeviceGuard guard(impl_->device, true);
impl_->stream = std::make_unique<CudaStream>();
impl_->available = true;
}
}
@@ -878,10 +895,10 @@ void RotationScaleMergeGPU::SetPartialsLayout(int n_obs, int n_frames,
namespace {
template <typename T>
void UploadChunk(CudaDevicePtr<T> &dst, int offset, int count, const T *v) {
void UploadChunk(CudaDevicePtr<T> &dst, int offset, int count, const T *v, cudaStream_t s) {
if (count > 0)
CudaCheck(cudaMemcpy(dst.get() + offset, v, size_t(count) * sizeof(T),
cudaMemcpyHostToDevice), "upload chunk");
CopyAndWait(dst.get() + offset, v, size_t(count) * sizeof(T), cudaMemcpyHostToDevice, s,
"upload chunk");
}
}
@@ -889,34 +906,34 @@ void RotationScaleMergeGPU::SetObsField(ObsField f, int offset, int count, const
DeviceGuard guard(impl_->device, impl_->available);
auto &d = *impl_;
switch (f) {
case ObsField::I: UploadChunk(d.I, offset, count, v); break;
case ObsField::Sigma: UploadChunk(d.sigma, offset, count, v); break;
case ObsField::PrescalingCorr: UploadChunk(d.prescaling_corr, offset, count, v); break;
case ObsField::Partiality: UploadChunk(d.partiality, offset, count, v); break;
case ObsField::Zeta: UploadChunk(d.zeta, offset, count, v); break;
case ObsField::Corr0: UploadChunk(d.corr, offset, count, v); break;
case ObsField::Bkg: UploadChunk(d.bkg, offset, count, v); break;
case ObsField::VarBkg: UploadChunk(d.var_bkg, offset, count, v); break;
case ObsField::ImageNumber: UploadChunk(d.image_number, offset, count, v); break;
case ObsField::D: UploadChunk(d.d_obs, offset, count, v); break;
case ObsField::Px: UploadChunk(d.px_obs, offset, count, v); break;
case ObsField::Py: UploadChunk(d.py_obs, offset, count, v); break;
case ObsField::I: UploadChunk(d.I, offset, count, v, impl_->s()); break;
case ObsField::Sigma: UploadChunk(d.sigma, offset, count, v, impl_->s()); break;
case ObsField::PrescalingCorr: UploadChunk(d.prescaling_corr, offset, count, v, impl_->s()); break;
case ObsField::Partiality: UploadChunk(d.partiality, offset, count, v, impl_->s()); break;
case ObsField::Zeta: UploadChunk(d.zeta, offset, count, v, impl_->s()); break;
case ObsField::Corr0: UploadChunk(d.corr, offset, count, v, impl_->s()); break;
case ObsField::Bkg: UploadChunk(d.bkg, offset, count, v, impl_->s()); break;
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;
}
}
void RotationScaleMergeGPU::SetObsFrame(int offset, int count, const int32_t *frame) {
DeviceGuard guard(impl_->device, impl_->available);
UploadChunk(impl_->frame, offset, count, frame);
UploadChunk(impl_->frame, offset, count, frame, impl_->s());
}
void RotationScaleMergeGPU::SetObsOnIce(int offset, int count, const uint8_t *on_ice) {
DeviceGuard guard(impl_->device, impl_->available);
UploadChunk(impl_->on_ice, offset, count, on_ice);
UploadChunk(impl_->on_ice, offset, count, on_ice, impl_->s());
}
void RotationScaleMergeGPU::SetObsClipped(int offset, int count, const uint8_t *clipped) {
DeviceGuard guard(impl_->device, impl_->available);
UploadChunk(impl_->clipped, offset, count, clipped);
UploadChunk(impl_->clipped, offset, count, clipped, impl_->s());
}
void RotationScaleMergeGPU::SetGroups(int n_groups, const int32_t *group, const int32_t *group_perm,
@@ -934,8 +951,8 @@ void RotationScaleMergeGPU::SetGroups(int n_groups, const int32_t *group, const
void RotationScaleMergeGPU::SetCorr(const float *corr) {
DeviceGuard guard(impl_->device, impl_->available);
CudaCheck(cudaMemcpy(impl_->corr.get(), corr, size_t(impl_->n_obs) * sizeof(float),
cudaMemcpyHostToDevice), "upload corr");
CopyAndWait(impl_->corr.get(), corr, size_t(impl_->n_obs) * sizeof(float),
cudaMemcpyHostToDevice, impl_->s(), "upload corr");
}
void RotationScaleMergeGPU::ScalePartials(int iters, double min_partiality, bool /*has_d_min*/) {
@@ -943,42 +960,42 @@ void RotationScaleMergeGPU::ScalePartials(int iters, double min_partiality, bool
auto &d = *impl_;
// Reset per call: the host keeps the G of a frame across calls (RunScalingLoop), so a frame this
// call did not fit must read as unfitted, not as fitted with the value of the call before.
CudaCheck(cudaMemset(d.scaled.get(), 0, size_t(d.n_frames) * sizeof(uint8_t)), "memset scaled");
CudaCheck(cudaMemset(d.g.get(), 0, size_t(d.n_frames) * sizeof(double)), "memset g"); // unscaled g unused
CudaCheck(cudaMemsetAsync(d.scaled.get(), 0, size_t(d.n_frames) * sizeof(uint8_t), impl_->s()), "memset scaled");
CudaCheck(cudaMemsetAsync(d.g.get(), 0, size_t(d.n_frames) * sizeof(double), impl_->s()), "memset g"); // unscaled g unused
const int obs_blocks = (d.n_obs + BLK - 1) / BLK;
const int upd_blocks = std::min(65535, obs_blocks);
const int grp_blocks = std::min(65535, (d.n_groups + BLK - 1) / BLK);
for (int it = 0; it < iters; ++it) {
ReduceGroupMeansKernel<<<grp_blocks, BLK>>>(d.n_groups, min_partiality,
ReduceGroupMeansKernel<<<grp_blocks, BLK, 0, impl_->s()>>>(d.n_groups, min_partiality,
d.group_perm.get(), d.group_start.get(), d.group_count.get(),
d.I.get(), d.sigma.get(), d.partiality.get(), d.corr.get(), d.group_mean.get());
CudaCheck(cudaGetLastError(), "ReduceGroupMeansKernel launch");
PrepScaleObsKernel<<<obs_blocks, BLK>>>(d.n_obs, min_partiality, d.group.get(), d.partiality.get(), d.prescaling_corr.get(),
PrepScaleObsKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(d.n_obs, min_partiality, d.group.get(), d.partiality.get(), d.prescaling_corr.get(),
d.zeta.get(), d.on_ice.get(), d.group_mean.get(), d.sigma.get(), d.inv_sigma.get(),
d.sco_coeff.get(), d.sco_ok.get());
CudaCheck(cudaGetLastError(), "PrepScaleObsKernel launch");
FitPerFrameGKernel<<<d.n_frames, BLK>>>(d.n_frames, d.frame_start.get(), d.frame_count.get(),
FitPerFrameGKernel<<<d.n_frames, BLK, 0, impl_->s()>>>(d.n_frames, d.frame_start.get(), d.frame_count.get(),
d.I.get(), d.inv_sigma.get(), d.sco_coeff.get(), d.sco_ok.get(), nullptr, d.g.get(), d.scaled.get());
CudaCheck(cudaGetLastError(), "FitPerFrameGKernel launch");
UpdateCorrKernel<<<upd_blocks, BLK>>>(d.n_obs, d.frame.get(), d.prescaling_corr.get(), d.partiality.get(),
UpdateCorrKernel<<<upd_blocks, BLK, 0, impl_->s()>>>(d.n_obs, d.frame.get(), d.prescaling_corr.get(), d.partiality.get(),
d.g.get(), d.scaled.get(), d.corr.get());
}
CudaCheck(cudaGetLastError(), "kernel launch");
CudaCheck(cudaDeviceSynchronize(), "scale sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "scale sync");
}
void RotationScaleMergeGPU::GetCorr(float *corr_out) const {
DeviceGuard guard(impl_->device, impl_->available);
CudaCheck(cudaMemcpy(corr_out, impl_->corr.get(), size_t(impl_->n_obs) * sizeof(float),
cudaMemcpyDeviceToHost), "download corr");
CopyAndWait(corr_out, impl_->corr.get(), size_t(impl_->n_obs) * sizeof(float),
cudaMemcpyDeviceToHost, impl_->s(), "download corr");
}
void RotationScaleMergeGPU::GetG(double *g_out, uint8_t *scaled_out) const {
DeviceGuard guard(impl_->device, impl_->available);
CudaCheck(cudaMemcpy(g_out, impl_->g.get(), size_t(impl_->n_frames) * sizeof(double),
cudaMemcpyDeviceToHost), "download g");
CudaCheck(cudaMemcpy(scaled_out, impl_->scaled.get(), size_t(impl_->n_frames) * sizeof(uint8_t),
cudaMemcpyDeviceToHost), "download scaled");
CopyAndWait(g_out, impl_->g.get(), size_t(impl_->n_frames) * sizeof(double),
cudaMemcpyDeviceToHost, impl_->s(), "download g");
CopyAndWait(scaled_out, impl_->scaled.get(), size_t(impl_->n_frames) * sizeof(uint8_t),
cudaMemcpyDeviceToHost, impl_->s(), "download scaled");
}
void RotationScaleMergeGPU::SetFrameCellOk(const uint8_t *frame_cell_ok) {
@@ -1030,19 +1047,19 @@ void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partialit
const int grp_blocks = std::min(65535, (ng + BLK - 1) / BLK);
const int obs_blocks = std::min(65535, (nf + BLK - 1) / BLK);
MergeEmStatsKernel<<<grp_blocks, BLK>>>(p);
MergeSamplesKernel<<<obs_blocks, BLK>>>(nf, p);
MergeEmStatsKernel<<<grp_blocks, BLK, 0, impl_->s()>>>(p);
MergeSamplesKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(nf, p);
CudaCheck(cudaGetLastError(), "merge em/samples launch");
CudaCheck(cudaDeviceSynchronize(), "merge em/samples sync");
CudaCheck(cudaMemcpy(em_mean_out, d.m_em_mean.get(), size_t(ng) * sizeof(double),
cudaMemcpyDeviceToHost), "dl em_mean");
CudaCheck(cudaMemcpy(cnt_out, d.m_cnt.get(), size_t(ng) * sizeof(int32_t),
cudaMemcpyDeviceToHost), "dl cnt");
CudaCheck(cudaStreamSynchronize(impl_->s()), "merge em/samples sync");
CopyAndWait(em_mean_out, d.m_em_mean.get(), size_t(ng) * sizeof(double),
cudaMemcpyDeviceToHost, impl_->s(), "dl em_mean");
CopyAndWait(cnt_out, d.m_cnt.get(), size_t(ng) * sizeof(int32_t),
cudaMemcpyDeviceToHost, impl_->s(), "dl cnt");
if (nf > 0) {
CudaCheck(cudaMemcpy(s2_out, d.m_s2.get(), size_t(nf) * sizeof(double), cudaMemcpyDeviceToHost), "dl s2");
CudaCheck(cudaMemcpy(I2_out, d.m_I2.get(), size_t(nf) * sizeof(double), cudaMemcpyDeviceToHost), "dl I2");
CudaCheck(cudaMemcpy(dev2_out, d.m_dev2.get(), size_t(nf) * sizeof(double), cudaMemcpyDeviceToHost), "dl dev2");
CudaCheck(cudaMemcpy(valid_out, d.m_valid.get(), size_t(nf) * sizeof(uint8_t), cudaMemcpyDeviceToHost), "dl valid");
CopyAndWait(s2_out, d.m_s2.get(), size_t(nf) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl s2");
CopyAndWait(I2_out, d.m_I2.get(), size_t(nf) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl I2");
CopyAndWait(dev2_out, d.m_dev2.get(), size_t(nf) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl dev2");
CopyAndWait(valid_out, d.m_valid.get(), size_t(nf) * sizeof(uint8_t), cudaMemcpyDeviceToHost, impl_->s(), "dl valid");
}
}
@@ -1091,10 +1108,10 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_
p.rejected_obs = d.m_rejected.get();
const int grp_blocks = std::min(65535, (ng + BLK - 1) / BLK);
MergeAccumKernel<<<grp_blocks, BLK>>>(p);
MergeAccumKernel<<<grp_blocks, BLK, 0, impl_->s()>>>(p);
CudaCheck(cudaGetLastError(), "merge accum launch");
CudaCheck(cudaDeviceSynchronize(), "merge accum sync");
CudaCheck(cudaMemcpy(rejected_obs, d.m_rejected.get(), size_t(d.n_fulls) * sizeof(uint8_t), cudaMemcpyDeviceToHost),
CudaCheck(cudaStreamSynchronize(impl_->s()), "merge accum sync");
CopyAndWait(rejected_obs, d.m_rejected.get(), size_t(d.n_fulls) * sizeof(uint8_t), cudaMemcpyDeviceToHost, impl_->s(),
"dl rejected_obs");
}
@@ -1105,7 +1122,7 @@ void RotationScaleMergeGPU::MergeAccumRange(int g0, int n, double *swI, double *
DeviceGuard guard(impl_->device, impl_->available);
auto &d = *impl_;
auto dl = [&](void *h, const auto &s) {
CudaCheck(cudaMemcpy(h, s.get() + g0, size_t(n) * sizeof(*s.get()), cudaMemcpyDeviceToHost),
CopyAndWait(h, s.get() + g0, size_t(n) * sizeof(*s.get()), cudaMemcpyDeviceToHost, impl_->s(),
"dl accum"); };
dl(swI, d.a_swI); dl(sw, d.a_sw); dl(swIh0, d.a_swIh0); dl(swIh1, d.a_swIh1);
dl(swh0, d.a_swh0); dl(swh1, d.a_swh1); dl(swh_typ0, d.a_swht0); dl(swh_typ1, d.a_swht1);
@@ -1139,17 +1156,17 @@ void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, d
p.rejected_obs = d.m_rejected.get();
const int grp_blocks = std::min(65535, (ng + BLK - 1) / BLK);
MergeRmeasKernel<<<grp_blocks, BLK>>>(p);
MergeRmeasKernel<<<grp_blocks, BLK, 0, impl_->s()>>>(p);
CudaCheck(cudaGetLastError(), "merge rmeas launch");
CudaCheck(cudaDeviceSynchronize(), "merge rmeas sync");
CudaCheck(cudaMemcpy(absdev, d.r_absdev.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl absdev");
CudaCheck(cudaMemcpy(sumI, d.r_sumI.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl sumI");
CudaCheck(cudaMemcpy(wabsdev, d.r_wabsdev.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl wabsdev");
CudaCheck(cudaMemcpy(wsumI, d.r_wsumI.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl wsumI");
CudaCheck(cudaMemcpy(sumv, d.r_sumv.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl sumv");
CudaCheck(cudaMemcpy(sumv2, d.r_sumv2.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl sumv2");
CudaCheck(cudaMemcpy(n, d.r_n.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl rn");
CudaCheck(cudaMemcpy(nusable, d.r_nusable.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl rnusable");
CudaCheck(cudaStreamSynchronize(impl_->s()), "merge rmeas sync");
CopyAndWait(absdev, d.r_absdev.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl absdev");
CopyAndWait(sumI, d.r_sumI.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl sumI");
CopyAndWait(wabsdev, d.r_wabsdev.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl wabsdev");
CopyAndWait(wsumI, d.r_wsumI.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl wsumI");
CopyAndWait(sumv, d.r_sumv.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl sumv");
CopyAndWait(sumv2, d.r_sumv2.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost, impl_->s(), "dl sumv2");
CopyAndWait(n, d.r_n.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost, impl_->s(), "dl rn");
CopyAndWait(nusable, d.r_nusable.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost, impl_->s(), "dl rnusable");
}
void RotationScaleMergeGPU::SmoothCorr(const uint8_t *apply, const double *ratio) {
@@ -1158,10 +1175,10 @@ void RotationScaleMergeGPU::SmoothCorr(const uint8_t *apply, const double *ratio
d.Upload(d.smooth_apply, apply, d.n_frames);
d.Upload(d.smooth_ratio, ratio, d.n_frames);
const int blocks = std::min(65535, (d.n_obs + BLK - 1) / BLK);
SmoothCorrKernel<<<blocks, BLK>>>(d.n_obs, d.frame.get(), d.smooth_apply.get(),
SmoothCorrKernel<<<blocks, BLK, 0, impl_->s()>>>(d.n_obs, d.frame.get(), d.smooth_apply.get(),
d.smooth_ratio.get(), d.corr.get());
CudaCheck(cudaGetLastError(), "smooth corr launch");
CudaCheck(cudaDeviceSynchronize(), "smooth corr sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "smooth corr sync");
}
void RotationScaleMergeGPU::SmoothFullsCorr(const uint8_t *apply, const double *ratio) {
@@ -1171,23 +1188,23 @@ void RotationScaleMergeGPU::SmoothFullsCorr(const uint8_t *apply, const double *
d.Upload(d.smooth_apply, apply, d.n_frames);
d.Upload(d.smooth_ratio, ratio, d.n_frames);
const int blocks = std::min(65535, (d.n_fulls + BLK - 1) / BLK);
SmoothCorrKernel<<<blocks, BLK>>>(d.n_fulls, d.f_frame.get(), d.smooth_apply.get(),
SmoothCorrKernel<<<blocks, BLK, 0, impl_->s()>>>(d.n_fulls, d.f_frame.get(), d.smooth_apply.get(),
d.smooth_ratio.get(), d.f_corr.get());
CudaCheck(cudaGetLastError(), "smooth fulls corr launch");
CudaCheck(cudaDeviceSynchronize(), "smooth fulls corr sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "smooth fulls corr sync");
}
int64_t RotationScaleMergeGPU::FilterCorrByZeta(double min_zeta) {
DeviceGuard guard(impl_->device, impl_->available);
auto &d = *impl_;
CudaDevicePtr<unsigned long long> dropped = d.Alloc<unsigned long long>(1);
CudaCheck(cudaMemset(dropped.get(), 0, sizeof(unsigned long long)), "zero zeta drop count");
CudaCheck(cudaMemsetAsync(dropped.get(), 0, sizeof(unsigned long long), impl_->s()), "zero zeta drop count");
const int blocks = std::min(65535, (d.n_obs + BLK - 1) / BLK);
FilterZetaKernel<<<blocks, BLK>>>(d.n_obs, min_zeta, d.zeta.get(), d.corr.get(), dropped.get());
FilterZetaKernel<<<blocks, BLK, 0, impl_->s()>>>(d.n_obs, min_zeta, d.zeta.get(), d.corr.get(), dropped.get());
CudaCheck(cudaGetLastError(), "zeta filter launch");
CudaCheck(cudaDeviceSynchronize(), "zeta filter sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "zeta filter sync");
unsigned long long n = 0;
CudaCheck(cudaMemcpy(&n, dropped.get(), sizeof(unsigned long long), cudaMemcpyDeviceToHost),
CopyAndWait(&n, dropped.get(), sizeof(unsigned long long), cudaMemcpyDeviceToHost, impl_->s(),
"dl zeta drop count");
return static_cast<int64_t>(n);
}
@@ -1197,9 +1214,9 @@ void RotationScaleMergeGPU::FilterCorrByFrame(const uint8_t *reject) {
auto &d = *impl_;
d.Upload(d.filter_reject, reject, d.n_frames);
const int blocks = std::min(65535, (d.n_obs + BLK - 1) / BLK);
FilterFrameKernel<<<blocks, BLK>>>(d.n_obs, d.frame.get(), d.filter_reject.get(), d.corr.get());
FilterFrameKernel<<<blocks, BLK, 0, impl_->s()>>>(d.n_obs, d.frame.get(), d.filter_reject.get(), d.corr.get());
CudaCheck(cudaGetLastError(), "frame filter launch");
CudaCheck(cudaDeviceSynchronize(), "frame filter sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "frame filter sync");
}
void RotationScaleMergeGPU::ComputePartialCC(double min_partiality, double *cc_out, int64_t *cc_n_out) {
@@ -1208,19 +1225,19 @@ void RotationScaleMergeGPU::ComputePartialCC(double min_partiality, double *cc_o
const int grp_blocks = std::min(65535, (d.n_groups + BLK - 1) / BLK);
// Post-smooth group means (reuse the scaling reduce; reads the resident, smoothed corr), then the
// per-frame CC over the resident partials. Only the tiny per-frame cc/cc_n come back to the host.
ReduceGroupMeansKernel<<<grp_blocks, BLK>>>(d.n_groups, min_partiality,
ReduceGroupMeansKernel<<<grp_blocks, BLK, 0, impl_->s()>>>(d.n_groups, min_partiality,
d.group_perm.get(), d.group_start.get(), d.group_count.get(),
d.I.get(), d.sigma.get(), d.partiality.get(), d.corr.get(), d.group_mean.get());
CudaCheck(cudaGetLastError(), "ReduceGroupMeansKernel launch");
PerFrameCCKernel<<<d.n_frames, BLK>>>(d.n_frames, min_partiality,
PerFrameCCKernel<<<d.n_frames, BLK, 0, impl_->s()>>>(d.n_frames, min_partiality,
d.frame_start.get(), d.frame_count.get(), d.I.get(), d.sigma.get(), d.partiality.get(),
d.corr.get(), d.on_ice.get(), d.group.get(), d.group_mean.get(), d.cc.get(), d.cc_n.get());
CudaCheck(cudaGetLastError(), "partial CC launch");
CudaCheck(cudaDeviceSynchronize(), "partial CC sync");
CudaCheck(cudaMemcpy(cc_out, d.cc.get(), size_t(d.n_frames) * sizeof(double),
cudaMemcpyDeviceToHost), "download cc");
CudaCheck(cudaMemcpy(cc_n_out, d.cc_n.get(), size_t(d.n_frames) * sizeof(int64_t),
cudaMemcpyDeviceToHost), "download cc_n");
CudaCheck(cudaStreamSynchronize(impl_->s()), "partial CC sync");
CopyAndWait(cc_out, d.cc.get(), size_t(d.n_frames) * sizeof(double),
cudaMemcpyDeviceToHost, impl_->s(), "download cc");
CopyAndWait(cc_n_out, d.cc_n.get(), size_t(d.n_frames) * sizeof(int64_t),
cudaMemcpyDeviceToHost, impl_->s(), "download cc_n");
}
void RotationScaleMergeGPU::SetRawRuns(int n_runs, int n_perm, const int32_t *perm,
@@ -1247,8 +1264,8 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti
float max_frame_gap) {
DeviceGuard guard(impl_->device, impl_->available);
auto &d = *impl_;
CudaCheck(cudaMemcpy(d.rr_group.get(), rawrun_group, size_t(d.n_runs) * sizeof(int32_t),
cudaMemcpyHostToDevice), "upload rr_group");
CopyAndWait(d.rr_group.get(), rawrun_group, size_t(d.n_runs) * sizeof(int32_t),
cudaMemcpyHostToDevice, impl_->s(), "upload rr_group");
CombineParams p{};
p.n_runs = d.n_runs;
@@ -1267,13 +1284,13 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti
const int blocks = std::min(65535, (d.n_runs + BLK - 1) / BLK);
// Count pass: how many fulls each run emits.
CombineKernel<false><<<blocks, BLK>>>(p);
CombineKernel<false><<<blocks, BLK, 0, impl_->s()>>>(p);
CudaCheck(cudaGetLastError(), "combine count launch");
// Exclusive prefix sum on the host (deterministic) -> per-run output offset + total fulls.
std::vector<int32_t> nevents(d.n_runs);
CudaCheck(cudaMemcpy(nevents.data(), d.rr_nevents.get(), size_t(d.n_runs) * sizeof(int32_t),
cudaMemcpyDeviceToHost), "download nevents");
CopyAndWait(nevents.data(), d.rr_nevents.get(), size_t(d.n_runs) * sizeof(int32_t),
cudaMemcpyDeviceToHost, impl_->s(), "download nevents");
std::vector<int32_t> offset(d.n_runs);
int64_t acc = 0;
for (int r = 0; r < d.n_runs; ++r) { offset[r] = static_cast<int32_t>(acc); acc += nevents[r]; }
@@ -1292,8 +1309,8 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti
d.f_rlp = d.Alloc<float>(nf); d.f_zeta = d.Alloc<float>(nf);
d.f_inv_sigma = d.Alloc<double>(nf);
d.f_sco_coeff = d.Alloc<float>(nf); d.f_sco_ok = d.Alloc<uint8_t>(nf);
CudaCheck(cudaMemcpy(d.rr_offset.get(), offset.data(), size_t(d.n_runs) * sizeof(int32_t),
cudaMemcpyHostToDevice), "upload offset");
CopyAndWait(d.rr_offset.get(), offset.data(), size_t(d.n_runs) * sizeof(int32_t),
cudaMemcpyHostToDevice, impl_->s(), "upload offset");
p.rr_offset = d.rr_offset.get();
p.f_h = d.f_h.get(); p.f_k = d.f_k.get(); p.f_l = d.f_l.get();
@@ -1303,10 +1320,10 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti
p.f_var_bkg = d.f_var_bkg.get(); p.f_var_per_I = d.f_var_per_I.get();
p.f_on_ice = d.f_on_ice.get(); p.f_clipped = d.f_clipped.get();
if (d.n_fulls > 0) {
CombineKernel<true><<<blocks, BLK>>>(p);
CombineKernel<true><<<blocks, BLK, 0, impl_->s()>>>(p);
CudaCheck(cudaGetLastError(), "combine emit launch");
}
CudaCheck(cudaDeviceSynchronize(), "combine sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "combine sync");
return d.n_fulls;
}
@@ -1318,7 +1335,7 @@ void RotationScaleMergeGPU::GetFulls(int32_t *h, int32_t *k, int32_t *l, float *
const size_t n = static_cast<size_t>(dd.n_fulls);
if (n == 0) return;
auto dl = [&](void *dst, const void *src, size_t bytes) {
CudaCheck(cudaMemcpy(dst, src, bytes, cudaMemcpyDeviceToHost), "download fulls");
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));
@@ -1334,8 +1351,8 @@ void RotationScaleMergeGPU::GetFullsKeys(int32_t *frame, int32_t *group) const {
const auto &d = *impl_;
if (d.n_fulls == 0) return;
const size_t bytes = size_t(d.n_fulls) * sizeof(int32_t);
CudaCheck(cudaMemcpy(frame, d.f_frame.get(), bytes, cudaMemcpyDeviceToHost), "download f_frame");
CudaCheck(cudaMemcpy(group, d.f_group.get(), bytes, cudaMemcpyDeviceToHost), "download f_group");
CopyAndWait(frame, d.f_frame.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_frame");
CopyAndWait(group, d.f_group.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_group");
}
void RotationScaleMergeGPU::SetFullsFrameCSR(const int32_t *frame_perm, int n_perm,
@@ -1363,15 +1380,15 @@ void RotationScaleMergeGPU::ResetFullsScale() {
if (nf == 0) return;
const int obs_blocks = std::min(65535, (nf + BLK - 1) / BLK);
// Unity model: partiality/prescaling_corr/zeta = 1 so coeff = mean; corr starts at 1.
FillKernel<<<obs_blocks, BLK>>>(d.f_corr.get(), nf, 1.0f);
FillKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(d.f_corr.get(), nf, 1.0f);
CudaCheck(cudaGetLastError(), "FillKernel launch");
FillKernel<<<obs_blocks, BLK>>>(d.f_partiality.get(), nf, 1.0f);
FillKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(d.f_partiality.get(), nf, 1.0f);
CudaCheck(cudaGetLastError(), "FillKernel launch");
FillKernel<<<obs_blocks, BLK>>>(d.f_rlp.get(), nf, 1.0f);
FillKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(d.f_rlp.get(), nf, 1.0f);
CudaCheck(cudaGetLastError(), "FillKernel launch");
FillKernel<<<obs_blocks, BLK>>>(d.f_zeta.get(), nf, 1.0f);
FillKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(d.f_zeta.get(), nf, 1.0f);
CudaCheck(cudaGetLastError(), "FillKernel launch");
CudaCheck(cudaDeviceSynchronize(), "reset fulls scale sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "reset fulls scale sync");
}
void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality) {
@@ -1383,38 +1400,38 @@ void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality) {
const int grp_blocks = std::min(65535, (d.n_groups + BLK - 1) / BLK);
// Reset per call, as ScalePartials: the host keeps the G of a frame across calls.
CudaCheck(cudaMemset(d.scaled.get(), 0, size_t(d.n_frames) * sizeof(uint8_t)), "memset f scaled");
CudaCheck(cudaMemset(d.g.get(), 0, size_t(d.n_frames) * sizeof(double)), "memset f g");
CudaCheck(cudaMemsetAsync(d.scaled.get(), 0, size_t(d.n_frames) * sizeof(uint8_t), impl_->s()), "memset f scaled");
CudaCheck(cudaMemsetAsync(d.g.get(), 0, size_t(d.n_frames) * sizeof(double), impl_->s()), "memset f g");
for (int it = 0; it < iters; ++it) {
ReduceGroupMeansKernel<<<grp_blocks, BLK>>>(d.n_groups, min_partiality,
ReduceGroupMeansKernel<<<grp_blocks, BLK, 0, impl_->s()>>>(d.n_groups, min_partiality,
d.f_gperm.get(), d.f_gstart.get(), d.f_gcount.get(),
d.f_I.get(), d.f_sigma.get(), d.f_partiality.get(), d.f_corr.get(), d.group_mean.get());
CudaCheck(cudaGetLastError(), "ReduceGroupMeansKernel launch");
// Not grid-stride, so its grid has to cover every full - unlike the grid-stride kernels
// below, which the 65535 cap is there for. Capped, it would silently leave the tail of
// sco_coeff/sco_ok stale above 16.8M fulls.
PrepScaleObsKernel<<<(nf + BLK - 1) / BLK, BLK>>>(nf, min_partiality, d.f_group.get(), d.f_partiality.get(),
PrepScaleObsKernel<<<(nf + BLK - 1) / BLK, BLK, 0, impl_->s()>>>(nf, min_partiality, d.f_group.get(), d.f_partiality.get(),
d.f_rlp.get(), d.f_zeta.get(), d.f_on_ice.get(), d.group_mean.get(),
d.f_sigma.get(), d.f_inv_sigma.get(), d.f_sco_coeff.get(), d.f_sco_ok.get());
CudaCheck(cudaGetLastError(), "PrepScaleObsKernel launch");
FitPerFrameGKernel<<<d.n_frames, BLK>>>(d.n_frames,
FitPerFrameGKernel<<<d.n_frames, BLK, 0, impl_->s()>>>(d.n_frames,
d.f_frame_start.get(), d.f_frame_count.get(), d.f_I.get(), d.f_inv_sigma.get(),
d.f_sco_coeff.get(), d.f_sco_ok.get(), d.f_frame_perm.get(), d.g.get(), d.scaled.get());
CudaCheck(cudaGetLastError(), "FitPerFrameGKernel launch");
UpdateCorrKernel<<<obs_blocks, BLK>>>(nf, d.f_frame.get(), d.f_rlp.get(), d.f_partiality.get(),
UpdateCorrKernel<<<obs_blocks, BLK, 0, impl_->s()>>>(nf, d.f_frame.get(), d.f_rlp.get(), d.f_partiality.get(),
d.g.get(), d.scaled.get(), d.f_corr.get());
}
CudaCheck(cudaGetLastError(), "scale fulls launch");
CudaCheck(cudaDeviceSynchronize(), "scale fulls sync");
CudaCheck(cudaStreamSynchronize(impl_->s()), "scale fulls sync");
}
void RotationScaleMergeGPU::GetFullsCorr(float *corr) const {
DeviceGuard guard(impl_->device, impl_->available);
const auto &d = *impl_;
if (d.n_fulls == 0) return;
CudaCheck(cudaMemcpy(corr, d.f_corr.get(), size_t(d.n_fulls) * sizeof(float),
cudaMemcpyDeviceToHost), "download f_corr");
CopyAndWait(corr, d.f_corr.get(), size_t(d.n_fulls) * sizeof(float),
cudaMemcpyDeviceToHost, impl_->s(), "download f_corr");
}
void RotationScaleMergeGPU::GetFullsPxPy(float *px, float *py) const {
@@ -1422,8 +1439,8 @@ void RotationScaleMergeGPU::GetFullsPxPy(float *px, float *py) const {
const auto &d = *impl_;
if (d.n_fulls == 0) return;
const size_t bytes = size_t(d.n_fulls) * sizeof(float);
CudaCheck(cudaMemcpy(px, d.f_px.get(), bytes, cudaMemcpyDeviceToHost), "download f_px");
CudaCheck(cudaMemcpy(py, d.f_py.get(), bytes, cudaMemcpyDeviceToHost), "download f_py");
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");
}
void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I) const {
@@ -1431,8 +1448,8 @@ void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I) c
const auto &d = *impl_;
if (d.n_fulls == 0) return;
const size_t bytes = size_t(d.n_fulls) * sizeof(float);
CudaCheck(cudaMemcpy(var_bkg, d.f_var_bkg.get(), bytes, cudaMemcpyDeviceToHost), "download f_var_bkg");
CudaCheck(cudaMemcpy(var_per_I, d.f_var_per_I.get(), bytes, cudaMemcpyDeviceToHost),
CopyAndWait(var_bkg, d.f_var_bkg.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_var_bkg");
CopyAndWait(var_per_I, d.f_var_per_I.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(),
"download f_var_per_I");
}
@@ -1440,6 +1457,6 @@ void RotationScaleMergeGPU::SetFullsCorr(const float *corr) {
DeviceGuard guard(impl_->device, impl_->available);
auto &d = *impl_;
if (d.n_fulls == 0) return;
CudaCheck(cudaMemcpy(d.f_corr.get(), corr, size_t(d.n_fulls) * sizeof(float),
cudaMemcpyHostToDevice), "upload f_corr");
CopyAndWait(d.f_corr.get(), corr, size_t(d.n_fulls) * sizeof(float),
cudaMemcpyHostToDevice, impl_->s(), "upload f_corr");
}
rugnux/Rugnux.cpp
+130 -17
View File
@@ -6637,6 +6637,28 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
const auto &rot_ss = experiment_.GetScalingSettings();
const bool is_rotation = experiment_.IsRotationIndexing(); // rotation indexing -> rotation scaling/merge
std::optional<RotationScaleMerge> rsm;
// How many times rsm has been ingested on this pass: a re-ingest follows a reindex or a cell
// change, after which a merge made from the first ingest is a merge of other indices.
int rsm_ingests = 0;
// Two P1 merges of this pass's integration made beside the rest of the tail, on an engine of
// their own: the all-observation arm of the space-group search (see there) and the P1
// cross-check (see where it is written). Neither reads anything the search, the in-symmetry
// merge or the analyses decide, so neither has to wait for them. The engine is ingested together
// with rsm, before any merge writes per-frame values back onto the outcomes, so the two start
// from the same state, and it writes nothing back itself. Used only while rsm was ingested once
// (rsm_ingest); otherwise both merges are made on rsm, as before.
struct P1MergesAhead {
DiffractionExperiment x;
Logger log = Logger::Buffered(); // what the engine logs
Logger all_observations_log = Logger::Buffered(); // ...up to the all-observation merge
std::optional<RotationScaleMerge> engine;
std::future<void> ingested;
std::future<RotationScaleMerge::Result> all_observations;
std::future<RotationScaleMerge::Result> crosscheck;
bool crosscheck_made = false;
int rsm_ingest = 0;
};
std::unique_ptr<P1MergesAhead> p1_ahead;
std::optional<PostRefineObservations> prepass_postrefine_obs;
// The rotation geometry post-refinement (see its call sites below for what it is for).
const auto post_refine_geometry = [&] {
@@ -6767,6 +6789,60 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Rotation scaling/merging (RotationScaleMerge) does not support "
"wedge refinement");
// The conditions of p1_crosscheck and of the all-observation arm below that are known here;
// the rest (a pass that turns out superseded, a search that finds no point group) only means
// a merge is made and not used.
const bool all_observations_ahead = !geometry_prepass && !experiment_.GetGemmiSpaceGroup().has_value()
&& rot_ss.GetSearchMinZeta() > 0.0;
const bool crosscheck_ahead = !geometry_prepass && write_files && config_.write_merged
&& config_.write_p1_crosscheck && result.consensus_cell
&& (!experiment_.GetGemmiSpaceGroup().has_value() || indexer->GetPredictionCentring() == 'P');
if ((all_observations_ahead || crosscheck_ahead) && config_.observation_dump_path.empty()) {
p1_ahead = std::make_unique<P1MergesAhead>();
p1_ahead->x = experiment_;
p1_ahead->rsm_ingest = 1;
p1_ahead->crosscheck_made = crosscheck_ahead;
std::promise<void> ingested;
std::promise<RotationScaleMerge::Result> all_observations;
p1_ahead->ingested = ingested.get_future();
if (all_observations_ahead)
p1_ahead->all_observations = all_observations.get_future();
p1_ahead->crosscheck = std::async(std::launch::async,
[&a = *p1_ahead, &outcomes = indexer->GetIntegrationOutcome(),
cell = result.consensus_cell, iter = static_cast<int>(config_.scaling_iter),
nthreads = config_.nthreads, ingested = std::move(ingested),
all_observations = std::move(all_observations), all_observations_ahead,
crosscheck_ahead]() mutable -> RotationScaleMerge::Result {
try {
a.engine.emplace(a.x, outcomes, cell, iter, nthreads, a.log);
a.engine->SetWriteBackPerFrameScale(false);
a.engine->Ingest();
} catch (...) {
ingested.set_exception(std::current_exception());
throw;
}
ingested.set_value();
// Ingested in the group rsm was, merged in P1 - as both merges on rsm are.
a.x.SpaceGroupNumber(1);
if (all_observations_ahead) {
try {
a.engine->SetSearchMinZeta(0.0);
auto merged = a.engine->Run(/*for_search=*/true, /*full_stats=*/true,
/*measure_cc_before_corrections=*/false);
a.all_observations_log = a.log;
a.log = Logger::Buffered();
all_observations.set_value(std::move(merged));
} catch (...) {
all_observations.set_exception(std::current_exception());
throw;
}
}
if (!crosscheck_ahead)
return {};
return a.engine->Run(/*for_search=*/false, /*full_stats=*/true,
/*measure_cc_before_corrections=*/false);
});
}
// A reference MTZ is allowed for rotation: it fixes the space group / cell (on the CLI) and
// resolves the indexing ambiguity (below), but is NOT used to scale - the rotation merge stays
// self-consistent.
@@ -6774,6 +6850,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
++rsm_ingests;
if (p1_ahead)
p1_ahead->ingested.get();
}
// The geometry pre-pass reads the per-image reflections exactly twice more: they were just
@@ -7108,9 +7187,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// all-observation arm on both twin gates and promoted to R32 by the filtered one. So the
// losing arm's refusal is logged and carried to the report below; the rule is unchanged.
if (rsm && rsm->GetSearchMinZeta() > 0.0 && !sg_search.point_group_hm.empty()) {
std::optional<RotationScaleMerge::Result> ahead;
if (p1_ahead && p1_ahead->all_observations.valid() && p1_ahead->rsm_ingest == rsm_ingests) {
ahead = p1_ahead->all_observations.get();
p1_ahead->all_observations_log.ReplayInto(logger);
}
const double zeta = rsm->GetSearchMinZeta();
rsm->SetSearchMinZeta(0.0);
auto sm_all = scale_and_merge("P1, all observations", true);
auto sm_all = scale_and_merge("P1, all observations", true, false, std::move(ahead));
rsm->SetSearchMinZeta(zeta);
merged_filtered_isa = sg_opts.merge_isa; // the filtered arm's, before it is replaced
sg_opts.merge_isa = result.error_model_isa; // this arm's own error model
@@ -7789,6 +7873,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
++rsm_ingests;
}
const auto &uc = *result.consensus_cell;
logger.Info("{} names a cell of {}x the volume of its reference setting: the "
@@ -7984,6 +8069,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
result.consensus_cell, static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
++rsm_ingests;
}
}
}
@@ -8011,6 +8097,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
++rsm_ingests;
}
}
experiment_.SetSpaceGroup(sg);
@@ -8344,6 +8431,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
++rsm_ingests;
phase("Re-merging in the reference's frame");
sm = scale_and_merge(moved->short_name(), false);
const auto &uc = *result.consensus_cell;
@@ -8389,6 +8477,27 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
const auto &twin_sg_opt = experiment_.GetGemmiSpaceGroup();
const gemmi::SpaceGroup *twin_sg = twin_sg_opt ? &*twin_sg_opt : nullptr;
// Diffraction anisotropy (see where it is reported, below), made beside the analyses that come
// before it: it reads the merge, the integrated observations and the per-frame scales the merge
// wrote back, none of which they change, and most of it is gathering the observations.
std::future<decltype(sm.statistics.anisotropy)> anisotropy;
if (!geometry_prepass && !superseded && result.consensus_cell) {
AnisotropyRunInfo aniso_run;
if (sm.statistics.sweep_quality.measured && sm.statistics.sweep_quality.sweep_deg > 0.0f)
aniso_run.observed_rotation_deg = sm.statistics.sweep_quality.sweep_deg;
aniso_run.dose_term_in_scale_model = experiment_.GetScalingSettings().GetCorrectionSurfaces();
aniso_run.radiation_damage_relative_b = sm.statistics.radiation_damage_delta_b;
const float wedge_deg = experiment_.GetGoniometer() ? experiment_.GetGoniometer()->GetWedge_deg() : 0.0f;
anisotropy = std::async(std::launch::async,
[&, aniso_run, wedge_deg, cell = *result.consensus_cell,
rotation = experiment_.IsRotationIndexing()] {
return AnalyzeAnisotropy(sm.merged,
ScaledObservations(indexer->GetIntegrationOutcome(), rotation, twin_sg,
wedge_deg, 0.5, config_.nthreads),
cell, twin_sg, aniso_run);
});
}
// Not on the geometry pre-pass, nor on a superseded one: the analysis goes into that pass's
// statistics text and its written reflections, and neither survives the run. The promotion flag
// below is a different thing - it is what the SEARCH did, the second pass reads it, and it is
@@ -8619,21 +8728,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// unmerged observations, because a merge has exact Laue symmetry by construction and the
// tensor directions the symmetry forbids - the only place a dataset measures its own
// systematic error - are identically zero in it.
if (result.consensus_cell) {
AnisotropyRunInfo aniso_run;
if (sm.statistics.sweep_quality.measured && sm.statistics.sweep_quality.sweep_deg > 0.0f)
aniso_run.observed_rotation_deg = sm.statistics.sweep_quality.sweep_deg;
aniso_run.dose_term_in_scale_model =
experiment_.GetScalingSettings().GetCorrectionSurfaces();
aniso_run.radiation_damage_relative_b = sm.statistics.radiation_damage_delta_b;
sm.statistics.anisotropy = AnalyzeAnisotropy(
sm.merged,
ScaledObservations(indexer->GetIntegrationOutcome(),
experiment_.IsRotationIndexing(), twin_sg,
experiment_.GetGoniometer()
? experiment_.GetGoniometer()->GetWedge_deg() : 0.0f,
0.5, config_.nthreads),
*result.consensus_cell, twin_sg, aniso_run);
if (anisotropy.valid()) {
sm.statistics.anisotropy = anisotropy.get();
stats_text << AnisotropyToText(sm.statistics.anisotropy) << "\n";
}
}
@@ -8831,6 +8927,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
&& !geometry_prepass && !superseded && config_.write_p1_crosscheck
&& p1_integration_complete && is_rotation;
std::optional<RotationScaleMerge::Result> p1_merged_early;
const bool p1_ahead_usable = p1_crosscheck && p1_ahead && p1_ahead->crosscheck_made
&& p1_ahead->rsm_ingest == rsm_ingests;
// Model validation runs BEFORE the reflection files are written, because it is what settles the
// frame they are written in: the enantiomorph, which merged intensities cannot choose, and -
@@ -8865,12 +8963,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// P1 merge afterwards, through merge_to_written, as it did when it was made below. The
// validation's log lines are held and printed as one block once it is done.
ModelValidationResult validation;
if (p1_crosscheck && rsm) {
if (p1_crosscheck && rsm && !p1_ahead_usable) {
Logger held = Logger::Buffered();
auto pending = std::async(std::launch::async, validate, std::ref(held));
experiment_.SpaceGroupNumber(1);
rsm->SetWriteBackPerFrameScale(false);
p1_merged_early = rsm->Run(/*for_search=*/false, /*full_stats=*/true,
/*measure_cc_before_corrections=*/false);
rsm->SetWriteBackPerFrameScale(true);
experiment_.SetSpaceGroup(data_sg);
validation = pending.get();
held.ReplayInto(logger);
@@ -9046,6 +9146,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
const double em_a = result.error_model_a;
const double em_b = result.error_model_b;
const auto res_fit = result.resolution_fit_A;
const int iter_partials = result.scaling_iterations_partials;
const int iter_fulls = result.scaling_iterations_fulls;
const bool converged = result.scaling_converged;
// The unmerged MTZ below does not depend on this merge, so it is built meanwhile, from
// the experiment as it stands in the determined group. The merge rewrites each image's
// mosaicity, which the file's batch headers carry, so that is filled in only after it.
@@ -9062,7 +9165,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// Both the merge and the MTZ read the group from the experiment, so it is set for
// the whole of it and restored after.
experiment_.SpaceGroupNumber(1);
if (p1_ahead_usable) {
p1_merged_early = p1_ahead->crosscheck.get();
p1_ahead->log.ReplayInto(logger);
}
if (!p1_merged_early)
rsm->SetWriteBackPerFrameScale(false);
auto p1 = scale_and_merge("P1 cross-check", false, false, std::move(p1_merged_early));
rsm->SetWriteBackPerFrameScale(true);
// The scaler still holds the observations in the indexing they were merged in, so every
// relabelling since (the written setting, the model's indexing) is applied to this merge
// too - or it would describe the dataset on other axes than the merged output beside it,
@@ -9135,6 +9245,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
result.error_model_a = em_a;
result.error_model_b = em_b;
result.resolution_fit_A = res_fit;
result.scaling_iterations_partials = iter_partials;
result.scaling_iterations_fulls = iter_fulls;
result.scaling_converged = converged;
if (determined != nullptr && determined->number > 1)
logger.Info("P1 cross-check dataset written to {} ({} unique reflections): the "
"same observations merged in P1 instead of {}, so a wrong space group "