From fb0272023e1e70df8e07e8cfc1048625475ad7fe Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 31 Jul 2026 14:53:17 +0200 Subject: [PATCH] scale_merge: apply outlier rejection to the anomalous split on the GPU path The GPU merge kernel rejects outliers on the device and keeps a per-full flag there, but only returned the per-group counts. The host array the CPU path fills stayed all zero, and the anomalous I(+)/I(-) accumulator is host-side and unconditional - so with --reject-outliers and a GPU present, the observations the merged IMEAN dropped were still averaged into I(+) and I(-). The same command on a CPU-only host excluded them: the exported anomalous differences depended on whether a GPU was there. R_meas was unaffected, having its own device-side path that reads the flags in place. MergeAccum now hands the per-full flags back so every host-side reduction sees the same rejections. The comment claiming reject_outliers was excluded from the GPU path was never true. Co-Authored-By: Claude Opus 5 (1M context) --- image_analysis/scale_merge/RotationScaleMerge.cpp | 10 ++++++---- image_analysis/scale_merge/RotationScaleMergeGPU.cu | 4 +++- image_analysis/scale_merge/RotationScaleMergeGPU.h | 6 ++++-- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index ffc7281e..37d8313d 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1456,7 +1456,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // The em-stats / samples / merge-accumulate / R_meas reductions run on the resident, scaled fulls // (their group CSR is still on the device from scale-fulls) when fulls_resident; the host keeps the - // I2-sort, the (a,b) fit, the export and the statistics. reject_outliers is excluded upstream. + // I2-sort, the (a,b) fit, the export and the statistics. Outlier rejection runs on the device too, + // and MergeAccum hands the per-full flags back so the host-side reductions see the same rejections. bool use_gpu_merge = false; #ifdef JFJOCH_USE_CUDA use_gpu_merge = fulls_resident && !fulls.empty(); @@ -1656,7 +1657,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool struct Accum { double swI = 0, sw = 0, swIh[2] = {0, 0}, swh[2] = {0, 0}; size_t nh[2] = {0, 0}; float d = NAN; }; std::vector acc(n_groups); size_t reject_count = 0; - std::vector rejected_obs(fulls.size(), 0); // per-full outlier-rejected flag (mirrors the GPU path) + std::vector rejected_obs(fulls.size(), 0); // per-full outlier-rejected flag (both paths) bool did_gpu_acc = false; #ifdef JFJOCH_USE_CUDA if (use_gpu_merge) { @@ -1666,7 +1667,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool gpu_->MergeAccum(error_model_a, error_model_b, error_model_active, reject_outliers, reject_nsigma, reject_median.data(), aswI.data(), asw.data(), aswIh0.data(), aswIh1.data(), - aswh0.data(), aswh1.data(), anh0.data(), anh1.data(), ad.data(), arej.data()); + aswh0.data(), aswh1.data(), anh0.data(), anh1.data(), ad.data(), arej.data(), + rejected_obs.data()); for (int g = 0; g < n_groups; ++g) { Accum &a = acc[g]; a.swI = aswI[g]; a.sw = asw[g]; a.swIh[0] = aswIh0[g]; a.swIh[1] = aswIh1[g]; @@ -1857,7 +1859,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool anom.reserve(result.merged.size() * 2 + 1); for (const auto &o : fulls) { if (!usable_merge(o)) continue; - if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge (CPU path) + if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge const HKLKey ak = anom_keygen(o.h, o.k, o.l); const int hand = ak.plus ? 0 : 1; const float I_corr = o.I * o.corr; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 9b8123f4..676f3243 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -788,7 +788,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ bool reject_outliers, double reject_nsigma, const float *reject_median, double *swI, double *sw, double *swIh0, double *swIh1, double *swh0, double *swh1, int32_t *nh0, int32_t *nh1, double *d_out, - int32_t *rejected) { + int32_t *rejected, uint8_t *rejected_obs) { auto &d = *impl_; const int ng = d.n_groups; d.a_swI = CudaDevicePtr(std::max(1, ng)); d.a_sw = CudaDevicePtr(std::max(1, ng)); @@ -829,6 +829,8 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ CudaCheck(cudaMemcpy(nh0, d.a_nh0.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl nh0"); CudaCheck(cudaMemcpy(nh1, d.a_nh1.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl nh1"); CudaCheck(cudaMemcpy(rejected, d.a_rejected.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl rej"); + CudaCheck(cudaMemcpy(rejected_obs, d.m_rejected.get(), size_t(d.n_fulls) * sizeof(uint8_t), cudaMemcpyDeviceToHost), + "dl rejected_obs"); } void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, double *sumI, diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 79a0bcaa..b5c3d37a 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -73,12 +73,14 @@ public: // Per-group merge accumulators (inv-var sums + deterministic half-sets, error-model-corrected sigma // from a/b). Outputs length n_groups; rejected[g] counts outliers dropped (reject_median uploaded, NAN - // where none). Requires MergeEmSamples first (em_mean resident). + // where none). rejected_obs is the per-full flag (length n_fulls): the host needs it for the reductions + // it still does itself, above all the anomalous I(+)/I(-) split. + // Requires MergeEmSamples first (em_mean resident). void MergeAccum(double error_model_a, double error_model_b, bool error_model_active, bool reject_outliers, double reject_nsigma, const float *reject_median, double *swI, double *sw, double *swIh0, double *swIh1, double *swh0, double *swh1, int32_t *nh0, int32_t *nh1, double *d_out, - int32_t *rejected); + int32_t *rejected, uint8_t *rejected_obs); // Per-group R_meas accumulators (sum|I_corr-merged_I|, sum_I, n, and the usable count for per-shell // total_observations); merged_I is uploaded. All arrays length n_groups.