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

scale_merge: apply outlier rejection to the anomalous split on the GPU path

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-31 14:53:17 +02:00
co-authored by Claude Opus 5
parent f8beee7e87
commit fb0272023e
3 changed files with 13 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+6 -4
View File
@@ -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<Accum> acc(n_groups);
size_t reject_count = 0;
std::vector<char> rejected_obs(fulls.size(), 0); // per-full outlier-rejected flag (mirrors the GPU path)
std::vector<uint8_t> 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;
image_analysis/scale_merge/RotationScaleMergeGPU.cu
+3 -1
View File
@@ -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<double>(std::max(1, ng)); d.a_sw = CudaDevicePtr<double>(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,
image_analysis/scale_merge/RotationScaleMergeGPU.h
+4 -2
View File
@@ -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.