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.