diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index c32a166f0..6eff70b3e 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -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 diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index ea1947c8f..79278560f 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -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 // >= 1 without cutting the final in-symmetry merge. Reset to the manual limit afterwards. void SetDMinLimit(std::optional d_min_A) { d_min_limit = d_min_A; } @@ -413,6 +419,7 @@ private: std::unique_ptr 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 diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 8ed922149..6e54fb380 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -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 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 &dst, const T *src, int n) const { dst = Alloc(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()) // 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(); impl_->available = true; } } @@ -878,10 +895,10 @@ void RotationScaleMergeGPU::SetPartialsLayout(int n_obs, int n_frames, namespace { template - void UploadChunk(CudaDevicePtr &dst, int offset, int count, const T *v) { + void UploadChunk(CudaDevicePtr &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<<>>(d.n_groups, min_partiality, + ReduceGroupMeansKernel<<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<<>>(d.n_obs, min_partiality, d.group.get(), d.partiality.get(), d.prescaling_corr.get(), + PrepScaleObsKernel<<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, d.frame_start.get(), d.frame_count.get(), + FitPerFrameGKernel<<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<<>>(d.n_obs, d.frame.get(), d.prescaling_corr.get(), d.partiality.get(), + UpdateCorrKernel<<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<<>>(p); - MergeSamplesKernel<<>>(nf, p); + MergeEmStatsKernel<<s()>>>(p); + MergeSamplesKernel<<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<<>>(p); + MergeAccumKernel<<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<<>>(p); + MergeRmeasKernel<<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<<>>(d.n_obs, d.frame.get(), d.smooth_apply.get(), + SmoothCorrKernel<<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<<>>(d.n_fulls, d.f_frame.get(), d.smooth_apply.get(), + SmoothCorrKernel<<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 dropped = d.Alloc(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<<>>(d.n_obs, min_zeta, d.zeta.get(), d.corr.get(), dropped.get()); + FilterZetaKernel<<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(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<<>>(d.n_obs, d.frame.get(), d.filter_reject.get(), d.corr.get()); + FilterFrameKernel<<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<<>>(d.n_groups, min_partiality, + ReduceGroupMeansKernel<<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, min_partiality, + PerFrameCCKernel<<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<<>>(p); + CombineKernel<<s()>>>(p); CudaCheck(cudaGetLastError(), "combine count launch"); // Exclusive prefix sum on the host (deterministic) -> per-run output offset + total fulls. std::vector nevents(d.n_runs); - CudaCheck(cudaMemcpy(nevents.data(), d.rr_nevents.get(), size_t(d.n_runs) * sizeof(int32_t), - cudaMemcpyDeviceToHost), "download nevents"); + CopyAndWait(nevents.data(), d.rr_nevents.get(), size_t(d.n_runs) * sizeof(int32_t), + cudaMemcpyDeviceToHost, impl_->s(), "download nevents"); std::vector offset(d.n_runs); int64_t acc = 0; for (int r = 0; r < d.n_runs; ++r) { offset[r] = static_cast(acc); acc += nevents[r]; } @@ -1292,8 +1309,8 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti d.f_rlp = d.Alloc(nf); d.f_zeta = d.Alloc(nf); d.f_inv_sigma = d.Alloc(nf); d.f_sco_coeff = d.Alloc(nf); d.f_sco_ok = d.Alloc(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<<>>(p); + CombineKernel<<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(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<<>>(d.f_corr.get(), nf, 1.0f); + FillKernel<<s()>>>(d.f_corr.get(), nf, 1.0f); CudaCheck(cudaGetLastError(), "FillKernel launch"); - FillKernel<<>>(d.f_partiality.get(), nf, 1.0f); + FillKernel<<s()>>>(d.f_partiality.get(), nf, 1.0f); CudaCheck(cudaGetLastError(), "FillKernel launch"); - FillKernel<<>>(d.f_rlp.get(), nf, 1.0f); + FillKernel<<s()>>>(d.f_rlp.get(), nf, 1.0f); CudaCheck(cudaGetLastError(), "FillKernel launch"); - FillKernel<<>>(d.f_zeta.get(), nf, 1.0f); + FillKernel<<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<<>>(d.n_groups, min_partiality, + ReduceGroupMeansKernel<<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, + FitPerFrameGKernel<<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<<>>(nf, d.f_frame.get(), d.f_rlp.get(), d.f_partiality.get(), + UpdateCorrKernel<<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"); } diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index ca1dab038..3ef633636 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -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 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 engine; + std::future ingested; + std::future all_observations; + std::future crosscheck; + bool crosscheck_made = false; + int rsm_ingest = 0; + }; + std::unique_ptr p1_ahead; std::optional 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(); + p1_ahead->x = experiment_; + p1_ahead->rsm_ingest = 1; + p1_ahead->crosscheck_made = crosscheck_ahead; + std::promise ingested; + std::promise 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(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(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 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(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(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(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(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 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 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 "