diff --git a/common/Reflection.h b/common/Reflection.h index 3b74a695..1bba5180 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -22,6 +22,7 @@ struct Reflection { float d; float I; float bkg; + float var_bkg; // non-signal (background) part of sigma^2, carried to the merge float sigma; float dist_ewald; float rlp; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index a4eaf5b5..c1e3b629 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -53,7 +53,6 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe r2_sq = r2 * r2; r3 = settings.GetR3(); r3_sq = r3 * r3; - min_sigma_ratio = settings.GetMinimumSigmaInRegardsToI(); R = static_cast(std::ceil(r2)); G = 2 * R + 1; GG = G * G; @@ -131,6 +130,7 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, struct Rough { double I = 0.0, sigma = NAN, bkg = 0.0, obs_x = 0.0, obs_y = 0.0; double bkg_var = 0.0; // variance of the background ESTIMATE itself, bkg / n_bkg + double var_bkg = 0.0; // total non-signal variance carried to the merge int64_t I_sum = 0; // kept so I can be rebuilt after the radial background correction int n_inner = 0; int r_bin = 0; // rounded distance from the beam centre, indexes the radial curve @@ -194,10 +195,12 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset. out.bkg_var = out.bkg / n_bkg_used; + out.var_bkg = static_cast(n_inner) * out.bkg + + static_cast(n_inner) * n_inner * out.bkg_var; out.I_sum = I_sum; out.n_inner = static_cast(n_inner); const double var_bkg_term = static_cast(n_inner) * n_inner * out.bkg_var; - out.sigma = std::max(1.0, out.I * min_sigma_ratio); + out.sigma = 1.0; if (I_sum > 0) { out.sigma = std::max(out.sigma, std::sqrt(static_cast(I_sum) + var_bkg_term)); out.obs_x = static_cast(I_sum_x) / static_cast(I_sum); @@ -244,7 +247,8 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const auto &rh = rough[i]; if (!rh.ok) continue; results[i] = {static_cast(rh.I), static_cast(rh.sigma), static_cast(rh.bkg), - static_cast(rh.obs_x), static_cast(rh.obs_y), true, rh.has_obs}; + static_cast(rh.obs_x), static_cast(rh.obs_y), + static_cast(rh.var_bkg), true, rh.has_obs}; } return Finalize(predicted, npredicted, results, image_number); } @@ -409,15 +413,19 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // I = sum(P*(px-bkg)/v) / sum(P^2/v), so dI/dbkg = -wsum/den and the background estimate's // own error adds (wsum/den)^2 * var(bkg) - the same term the box sum was missing. double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var); + double var_bkg = std::max(0.0, 1.0 / den - std::max(0.0, I) + + (wsum / den) * (wsum / den) * rh.bkg_var); if (std::abs(I - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) { I = rh.I; sigma = rh.sigma; + var_bkg = rh.var_bkg; } // Carry the Pass-A box-sum intensity-weighted centroid (observed spot position) through the // profile path too - post-refinement uses it as the observed position (beam-centre / distance). results[i] = {static_cast(I), static_cast(sigma), static_cast(rh.bkg), - static_cast(rh.obs_x), static_cast(rh.obs_y), true, rh.has_obs}; + static_cast(rh.obs_x), static_cast(rh.obs_y), + static_cast(var_bkg), true, rh.has_obs}; } return Finalize(predicted, npredicted, results, image_number); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 9e01139d..c6046ac2 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -16,7 +16,6 @@ inline void cuda_err(cudaError_t val) { struct BraggGpuParams { int W, H; float r1_sq, r2, r2_sq, r3, r3_sq; - float min_sigma_ratio; int R, G, GG; float bkg_clip_nsigma; // high-side background sigma-clip multiplier (0 = no clip) int empirical; // ProfileEmpirical vs ProfileGaussian @@ -56,7 +55,7 @@ __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, B __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int32_t *img, const uint8_t *mask, BraggGpuParams p, int n, int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o, - float *bkgvar_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o, + float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o, uint8_t *hasobs_o, unsigned long long *invd2mm, float *isum_o, int *ninner_o, int *rbin_o, float *rad_sum, int *rad_cnt, int n_rad) { @@ -235,7 +234,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, // var(I) = Isum + n_inner^2 * bkg/n_bkg_used. Both engines must agree. const double bkg_var = bkg / (double) n_bkg_used; const double var_bkg_term = (double) s_ninner * (double) s_ninner * bkg_var; - double sigma = fmax(1.0, I * (double) p.min_sigma_ratio); + double sigma = 1.0; uint8_t hasobs = 0; double ox = 0.0, oy = 0.0; if (Isum > 0) { sigma = fmax(sigma, sqrt((double) Isum + var_bkg_term)); @@ -247,6 +246,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, cy_o[i] = (int) lroundf(cy); I_o[i] = (float) I; sigma_o[i] = (float) sigma; bkg_o[i] = (float) bkg; bkgvar_o[i] = (float) bkg_var; + varbkg_o[i] = (float) ((double) s_ninner * bkg + var_bkg_term); isum_o[i] = (float) Isum; ninner_o[i] = s_ninner; rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0); @@ -383,10 +383,10 @@ __global__ void radial_correct(const float *rad_sum, const int *rad_cnt, int n_r __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_seed, const float *sigma_seed, const float *bkg_a, - const float *bkgvar_a, const uint8_t *ok_a, + const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, const float *shell_P, const float *global_P, const float *shell_sigma2, const float *global_sigma2, const int *shell_n, - float *I_o, float *sigma_o, uint8_t *ok_o, BraggGpuParams p, int n) { + float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n) { const int i = blockIdx.x; if (i >= n) return; extern __shared__ float Pbuf[]; @@ -471,13 +471,15 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, // estimate's own error (see the CPU engine). const float wr = s_wsum / s_den; float I = s_I, sigma = sqrtf(1.0f / s_den + wr * wr * bkgvar_a[i]); + float var_bkg = fmaxf(0.0f, 1.0f / s_den - fmaxf(0.0f, I) + wr * wr * bkgvar_a[i]); // Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed // when the profile result diverges from it. if (fabsf(I - I_seed[i]) > (float) PROFILE_SUMMATION_MAX_NSIGMA * sigma_seed[i]) { I = I_seed[i]; sigma = sigma_seed[i]; + var_bkg = varbkg_seed[i]; } - I_o[i] = I; sigma_o[i] = sigma; ok_o[i] = 1; + I_o[i] = I; sigma_o[i] = sigma; varbkg_o[i] = var_bkg; ok_o[i] = 1; } else ok_o[i] = 0; } } @@ -546,6 +548,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { d_sigma = CudaDevicePtr(new_capacity); d_bkg = CudaDevicePtr(new_capacity); d_bkg_var = CudaDevicePtr(new_capacity); + d_var_bkg = CudaDevicePtr(new_capacity); d_isum = CudaDevicePtr(new_capacity); d_ninner = CudaDevicePtr(new_capacity); d_rbin = CudaDevicePtr(new_capacity); @@ -557,6 +560,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { h_px_x.resize(new_capacity); h_px_y.resize(new_capacity); h_d.resize(new_capacity); h_I.resize(new_capacity); h_sigma.resize(new_capacity); h_bkg.resize(new_capacity); + h_var_bkg.resize(new_capacity); h_obs_x.resize(new_capacity); h_obs_y.resize(new_capacity); h_ok.resize(new_capacity); h_has_obs.resize(new_capacity); capacity = new_capacity; @@ -588,7 +592,6 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu BraggGpuParams p{ .W = static_cast(xpixel), .H = static_cast(ypixel), .r1_sq = r1_sq, .r2 = r2, .r2_sq = r2_sq, .r3 = r3, .r3_sq = r3_sq, - .min_sigma_ratio = min_sigma_ratio, .R = R, .G = G, .GG = GG, .bkg_clip_nsigma = mode != IntegratorMode::BoxSum ? bkg_clip_nsigma : 0.0f, .empirical = empirical ? 1 : 0, @@ -613,7 +616,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_invd2, GG); mark_mask<<>>(d_px_x, d_px_y, d_mask, p, n); boxsum<<>>(d_px_x, d_px_y, d_d, img, d_mask, p, n, - d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_obs_x, d_obs_y, + d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y, d_ok, d_strong, d_has_obs, d_invd2, d_isum, d_ninner, d_rbin, d_rad_sum, d_rad_cnt, rad_n); @@ -633,14 +636,15 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, p); fit<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, - d_I, d_sigma, d_bkg, d_bkg_var, d_ok, d_shell_P, d_global_P, + d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, d_shell_n, - d_I, d_sigma, d_ok, p, n); + d_I, d_sigma, d_var_bkg, d_ok, p, n); } cuda_err(cudaMemcpyAsync(h_I.data(), d_I, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(h_sigma.data(), d_sigma, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(h_bkg.data(), d_bkg, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(h_var_bkg.data(), d_var_bkg, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(h_ok.data(), d_ok, sizeof(uint8_t) * npredicted, cudaMemcpyDeviceToHost, *stream)); // Pass A always fills the box-sum centroid (observed spot position), so copy it back in all modes - // post-refinement uses it as the observed position (beam-centre / distance). @@ -654,6 +658,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu results[i].I = h_I[i]; results[i].sigma = h_sigma[i]; results[i].bkg = h_bkg[i]; + results[i].var_bkg = h_var_bkg[i]; results[i].ok = true; if (h_has_obs[i]) { results[i].observed_x = h_obs_x[i]; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index f83d03fb..a9d0e4eb 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -29,7 +29,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // --- per-reflection device arrays (grown by EnsureCapacity) --- CudaDevicePtr d_px_x, d_px_y, d_d; CudaDevicePtr d_cx, d_cy; - CudaDevicePtr d_I, d_sigma, d_bkg, d_bkg_var, d_obs_x, d_obs_y; + CudaDevicePtr d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y; CudaDevicePtr d_isum; // box-sum raw sum, for the radial correction CudaDevicePtr d_ninner, d_rbin; CudaDevicePtr d_ok, d_strong, d_has_obs; @@ -51,7 +51,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // --- host staging (copied back once per frame) --- std::vector h_px_x, h_px_y, h_d; - std::vector h_I, h_sigma, h_bkg, h_obs_x, h_obs_y; + std::vector h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y; std::vector h_ok, h_has_obs; void EnsureCapacity(size_t n); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 4a99ea95..0db26709 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -82,7 +82,7 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id continue; auto hkl = generator(r); auto hkl_key = hkl.pack(); - sigma_corr = CorrectedSigma(I_corr, sigma_corr, r.image_scale_corr, hkl_key); + sigma_corr = CorrectedSigma(I_corr, sigma_corr, r.image_scale_corr, r.var_bkg, hkl_key); // Robust outlier rejection: drop this observation if it sits more than // reject_nsigma error-model sigmas from the reflection's median. Needs the active @@ -121,6 +121,7 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id } float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, float image_scale_corr, + float var_bkg, uint64_t hkl_key) const { if (!error_model_active) return sigma_corr; @@ -144,8 +145,7 @@ float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, float image_ // observations this correction exists to stop being mistreated. double a_var = static_cast(sigma_corr) * sigma_corr; if (scaling_settings.GetExpectedVarianceMerge()) { - const double bkg_var = std::max(0.0, a_var - static_cast(image_scale_corr) - * std::max(0.0, static_cast(I_corr))); + const double bkg_var = static_cast(image_scale_corr) * image_scale_corr * var_bkg; const double base = bkg_var + static_cast(image_scale_corr) * std::max(0.0, I_for_b); if (base > 0.0) a_var = base; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 7adf10eb..464ac152 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -128,6 +128,7 @@ class MergeOnTheFly { // using the per-observation I_i instead would over-weight down-fluctuated points. std::unordered_map error_model_mean_I; [[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, float image_scale_corr, + float var_bkg, uint64_t hkl_key) const; // Optional per-observation outlier rejection: drop observations whose corrected diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 225d8cbe..d73ace1f 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -226,7 +226,7 @@ void RotationScaleMerge::Ingest() { Obs obs{}; obs.h = r.h; obs.k = r.k; obs.l = r.l; obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; obs.rlp = r.rlp; - obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; + obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg; obs.px = r.predicted_x; obs.py = r.predicted_y; obs.image_number = r.image_number; obs.frame = o; @@ -306,15 +306,17 @@ void RotationScaleMerge::Ingest() { px(n), py(n); std::vector onice(n); std::vector frm(n); + std::vector vbkg(n); for (int i = 0; i < n; ++i) { const auto &o = partials[i]; I[i] = o.I; sigma[i] = o.sigma; rlp[i] = o.rlp; part[i] = o.partiality; zeta[i] = o.zeta; onice[i] = o.on_ice; frm[i] = o.frame; corr[i] = o.corr; - bkg[i] = o.bkg; img[i] = o.image_number; dd[i] = o.d; px[i] = o.px; py[i] = o.py; + bkg[i] = o.bkg; vbkg[i] = o.var_bkg; img[i] = o.image_number; dd[i] = o.d; + px[i] = o.px; py[i] = o.py; } gpu_->SetPartials(n, n_frames, I.data(), sigma.data(), rlp.data(), part.data(), zeta.data(), onice.data(), frm.data(), corr.data(), frame_start.data(), frame_count.data()); - gpu_->SetCombineInputs(bkg.data(), img.data(), dd.data(), px.data(), py.data()); + gpu_->SetCombineInputs(bkg.data(), vbkg.data(), img.data(), dd.data(), px.data(), py.data()); gpu_->SetRawRuns(static_cast(rawrun_start.size()), static_cast(perm.size()), perm.data(), rawrun_start.data(), rawrun_count.data(), rawrun_h.data(), rawrun_k.data(), rawrun_l.data()); @@ -1457,9 +1459,11 @@ void RotationScaleMerge::Combine() { const double corr = r2.corr; const double I_corr = pooled_I(r2) * corr; const double sigma_corr = static_cast(r2.sigma) * corr; - // max(0, I): a down-fluctuated partial carries no Poisson signal to remove, and - // removing a negative one inflates the background part instead of leaving it alone. - const double bkg_var = sigma_corr * sigma_corr - corr * std::max(0.0, I_corr); + // The non-signal variance as the integrator measured it. It used to be + // back-derived here as sigma^2 - I, which assumes sigma^2 = I + N exactly - true + // for a box sum, never true for a profile fit, and false for anything whose + // sigma was floored. + const double bkg_var = corr * corr * static_cast(r2.var_bkg); double var = std::max(0.0, bkg_var) + corr * std::max(0.0, F); if (!(var > 0.0)) var = sigma_corr * sigma_corr; const double w = 1.0 / var; @@ -1812,12 +1816,24 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool auto &g = gs[o.group]; g.sum += I_corr; g.sum_sq += I_corr * I_corr; g.sum_var += sigma_corr * sigma_corr; ++g.n; } + // The counting term to subtract. The estimator below has always used the RAW reported + // sigma^2, but the (a, b) fit immediately above concluded the counting variance is + // a*sigma^2. Where a < 1 the raw subtraction overshoots: the per-group systematic is driven + // to zero, the median lands on the boundary, and 1/b_asy reports an impossible I/sigma + // (measured: one dataset reporting 64.6 where every other statistic supports ~16, and + // flipping between 10.9 and 64.6 on consecutive merges of the same data). Opt-in for now + // (env JFJOCH_ISA_ASY_A) so the default path stays bit-identical while the two are + // batteried against each other -- an earlier attempt at this subtraction was rejected for + // over-claiming on a > 1 data, which this change would reintroduce. + const double asy_counting_scale = + std::getenv("JFJOCH_ISA_ASY_A") != nullptr ? error_model_a : 1.0; + // Per-group counting-subtracted fractional systematic variance, paired with the group's I/sigma. std::vector> group_scatter; // (systematic b^2, I/sigma) for (const auto &g : gs) { if (g.n < 5) continue; const double mean = g.sum / g.n; - const double counting = g.sum_var / g.n; + const double counting = asy_counting_scale * (g.sum_var / g.n); if (mean <= 0.0 || counting <= 0.0) continue; const double variance = (g.sum_sq - g.sum * g.sum / g.n) / (g.n - 1); group_scatter.push_back({std::max(variance - counting, 0.0) / (mean * mean), @@ -1832,6 +1848,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool double b_asy = asymptote_above(20.0, 100); // tight threshold on data that supports it if (b_asy <= 0.0) b_asy = asymptote_above(10.0, 50); // relaxed for weak / damaged data if (b_asy > 0.0) error_model_b_asymptotic = b_asy; + } // Guard a degenerate low-multiplicity fit: with too few symmetry equivalents both the (a, b) fit and // the per-group scatter collapse toward zero, and 1/b then reports an impossibly high asymptotic diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 15677ee5..1eb878fe 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -80,7 +80,7 @@ private: // scale-fulls/merge. Flat (not nested per image); a POD so the arrays translate straight to CUDA. struct Obs { int32_t h, k, l; - float I, sigma, d, rlp, partiality, zeta, delta_phi, bkg; + float I, sigma, d, rlp, partiality, zeta, delta_phi, bkg, var_bkg; float px = NAN, py = NAN; // predicted detector position (for the absorption surface; CPU path only) float image_number; // fractional frame position (for 3D-combine contiguity) int32_t frame; // index of the outcome whose per-frame scale G applies to this obs diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index b9b85b60..7f46bbfa 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -259,7 +259,7 @@ namespace { struct CombineParams { int n_runs; double min_partiality, capture_uncertainty_coeff, min_captured_fraction; - const float *I, *sigma, *corr, *partiality, *bkg, *image_number, *d, *px, *py; + const float *I, *sigma, *corr, *partiality, *bkg, *var_bkg, *image_number, *d, *px, *py; const int32_t *frame; const uint8_t *on_ice; const int32_t *perm, *rr_start, *rr_count, *rr_h, *rr_k, *rr_l, *rr_group; @@ -302,6 +302,7 @@ namespace { cursor = ev_end + 1; // Pass A: pooled background = mean of the event members' finite backgrounds. + double pooled_bkg = 0.0; int n_pool = 0; for (int m = ev_start; m <= ev_end; ++m) { @@ -356,7 +357,8 @@ namespace { const double sigma_corr = double(p.sigma[i]) * corr; // max(0, I): a down-fluctuated partial carries no Poisson signal to remove, and // removing a negative one inflates the background part instead of leaving it alone. - const double bkg_var = sigma_corr * sigma_corr - corr * Dmax(0.0, I_corr); + // The integrator's own non-signal variance (see the host combine). + const double bkg_var = corr * corr * (double) p.var_bkg[i]; double var = Dmax(0.0, bkg_var) + corr * Dmax(0.0, F); if (!(var > 0.0)) var = sigma_corr * sigma_corr; const double w = 1.0 / var; @@ -618,7 +620,7 @@ struct RotationScaleMergeGPU::Impl { double merge_min_part = 0.0; // combine: extra per-obs inputs + the one-time raw-hkl run layout - CudaDevicePtr bkg, image_number, d_obs, px_obs, py_obs; + CudaDevicePtr bkg, var_bkg, image_number, d_obs, px_obs, py_obs; int n_runs = 0, n_perm = 0; CudaDevicePtr perm, rr_start, rr_count, rr_h, rr_k, rr_l, rr_group; CudaDevicePtr rr_nevents, rr_nusable, rr_offset; @@ -871,10 +873,12 @@ void RotationScaleMergeGPU::ComputePartialCC(double min_partiality, double *cc_o cudaMemcpyDeviceToHost), "download cc_n"); } -void RotationScaleMergeGPU::SetCombineInputs(const float *bkg, const float *image_number, const float *d, +void RotationScaleMergeGPU::SetCombineInputs(const float *bkg, const float *var_bkg, + const float *image_number, const float *d, const float *px, const float *py) { auto &dd = *impl_; Upload(dd.bkg, bkg, dd.n_obs); + Upload(dd.var_bkg, var_bkg, dd.n_obs); Upload(dd.image_number, image_number, dd.n_obs); Upload(dd.d_obs, d, dd.n_obs); Upload(dd.px_obs, px, dd.n_obs); @@ -911,7 +915,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.capture_uncertainty_coeff = capture_uncertainty_coeff; p.min_captured_fraction = min_captured_fraction; p.I = d.I.get(); p.sigma = d.sigma.get(); p.corr = d.corr.get(); p.partiality = d.partiality.get(); - p.bkg = d.bkg.get(); p.image_number = d.image_number.get(); p.d = d.d_obs.get(); + p.bkg = d.bkg.get(); p.var_bkg = d.var_bkg.get(); p.image_number = d.image_number.get(); p.d = d.d_obs.get(); p.px = d.px_obs.get(); p.py = d.py_obs.get(); p.frame = d.frame.get(); p.on_ice = d.on_ice.get(); p.perm = d.perm.get(); p.rr_start = d.rr_start.get(); p.rr_count = d.rr_count.get(); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 0f8795b9..9aea3739 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -96,7 +96,8 @@ public: // The per-obs fields the combine needs on top of the scaling inputs (image-local bkg, fractional // frame position for event contiguity, resolution, and the predicted detector position px/py carried // to the full for the absorption surface). Uploaded once, alongside SetPartials. - void SetCombineInputs(const float *bkg, const float *image_number, const float *d, + void SetCombineInputs(const float *bkg, const float *var_bkg, + const float *image_number, const float *d, const float *px, const float *py); // The one-time raw-hkl run layout (space-group-independent): the (raw h,k,l, image_number)-sorted diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index f119efb7..574dd6ff 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -150,6 +150,8 @@ bool ReadReflectionsFromGroup(HDF5Object &file, auto int_sum = file.ReadOptVector(image_group_name + "/int_sum"); auto int_err = file.ReadOptVector(image_group_name + "/int_err"); auto bkg = file.ReadOptVector(image_group_name + "/background_mean"); + // Written since the merge stopped back-deriving it; older _process.h5 do not carry it. + auto var_bkg = file.ReadOptVector(image_group_name + "/background_variance"); auto lp = file.ReadOptVector(image_group_name + "/lp"); auto partiality = file.ReadOptVector(image_group_name + "/partiality"); auto phi = file.ReadOptVector(image_group_name + "/delta_phi"); @@ -177,6 +179,12 @@ bool ReadReflectionsFromGroup(HDF5Object &file, if (zeta.size() > i) zeta_val = zeta[i]; + // A file written before the merge carried this cannot say what the non-signal variance was; + // 0 leaves the combine with the signal term alone, which is what it had before. + float var_bkg_val = 0.0f; + if (var_bkg.size() > i) + var_bkg_val = var_bkg[i]; + float image_scale_corr_val = 1.0f; // Default is 1.0, if we don't know any better if (image_scale_corr.size() > i) image_scale_corr_val = image_scale_corr[i]; @@ -202,6 +210,7 @@ bool ReadReflectionsFromGroup(HDF5Object &file, .d = d.at(i), .I = int_sum.at(i), .bkg = bkg.at(i), + .var_bkg = var_bkg_val, .sigma = int_err.at(i), .rlp = lp_val, .partiality = partiality_val, diff --git a/writer/HDF5DataFilePluginReflection.cpp b/writer/HDF5DataFilePluginReflection.cpp index 5a35c797..f858d6c9 100644 --- a/writer/HDF5DataFilePluginReflection.cpp +++ b/writer/HDF5DataFilePluginReflection.cpp @@ -16,7 +16,8 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ std::vector h, k, l; std::vector I, sigma, d, lp; - std::vector image, phi, pred_x, pred_y, obs_x, obs_y, bkg, partiality, zeta, scale_factor; + std::vector image, phi, pred_x, pred_y, obs_x, obs_y, bkg, var_bkg, partiality, zeta, + scale_factor; h.reserve(msg.reflections.size()); k.reserve(msg.reflections.size()); @@ -29,6 +30,7 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ obs_x.reserve(msg.reflections.size()); obs_y.reserve(msg.reflections.size()); bkg.reserve(msg.reflections.size()); + var_bkg.reserve(msg.reflections.size()); lp.reserve(msg.reflections.size()); partiality.reserve(msg.reflections.size()); image.reserve(msg.reflections.size()); @@ -49,6 +51,7 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ obs_x.emplace_back(refl.observed_x); obs_y.emplace_back(refl.observed_y); bkg.emplace_back(refl.bkg); + var_bkg.emplace_back(refl.var_bkg); lp.emplace_back(1.0/refl.rlp); partiality.emplace_back(refl.partiality); phi.emplace_back(refl.delta_phi_deg); @@ -71,6 +74,7 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ image_group.SaveVector("int_sum", I); image_group.SaveVector("int_err", sigma); image_group.SaveVector("background_mean", bkg); + image_group.SaveVector("background_variance", var_bkg); image_group.SaveVector("observed_frame", image); image_group.SaveVector("lp", lp); image_group.SaveVector("partiality", partiality);