diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 191546b6..9d1db4f6 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -4,6 +4,8 @@ This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. * Space-group search: a screw axis is now claimed from how decisively its predicted-absent reflections are weaker than the rest of their own axial row, instead of from a minimum count of them, so a screw survives a sweep that recorded few axial reflections and is refused on a row too weak to decide either way. +* Bragg integration: the profile fit's `background_variance` now takes the fitted intensity itself out of the fit variance instead of `max(0, I)`, so a reflection that fluctuated below zero no longer reports a background variance two to three times too small and is no longer weighted up for it. +* Scaling: the rotation merge weights each combined full by its variance rebuilt at the reflection's mean intensity rather than by the full's own sigma, as the stills merge already did. * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** * Bragg integration: the background ring's outer radius default changes from 10 px to **13 px**, which roughly doubles the pixels behind each background estimate; the signal disk is unchanged. * Scaling: when too few reflections are strong enough to constrain the error model's systematic term **b**, it is now held at zero and **ISa is reported as unmeasured** rather than as the spurious value the fit would otherwise extrapolate. diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 2916dd84..84c5bb36 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -511,7 +511,13 @@ 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) + // var_bkg is the NON-SIGNAL part of that variance, and 1/den is the fit's variance taken at the + // fitted intensity, so the signal part to remove is I itself - not max(0, I). Clamping it leaves + // a down-fluctuated reflection's own (deflated) variance standing as its background variance, + // which is 2-3x too small; the merge then weights exactly the down-fluctuated observations up. + // The whole point of a separate var_bkg is that it does not move with the observation's own + // fluctuation, and 1/den - I is what does not (1/den grows with I one for one). + double var_bkg = std::max(0.0, 1.0 / den - I + (wsum / den) * (wsum / den) * rh.bkg_var); // The seed is a sum over the disk the box sum actually read, so when Exclude has taken pixels // out of both, the fit's full-profile intensity has to be scaled down to that same disk diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 49229085..986c272b 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -574,7 +574,8 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x // 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]); + // The signal part to remove is I itself, not max(0, I) - see the CPU engine. + float var_bkg = fmaxf(0.0f, 1.0f / s_den - 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. Exclude has taken pixels out of both, so the // full-profile intensity is scaled back to the disk the seed read; nothing excluded gives diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index f9e43a73..94439bcd 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1570,26 +1570,31 @@ void RotationScaleMerge::Combine() { } double F = sum_wI / sum_w; + // The full's variance as a function of intensity, var(I) = var_bkg_full + slope*I, so the + // merge can rebuild it at the reflection's mean instead of at this event's own answer. + double sum_wb = 0.0, sum_cwb = 0.0; for (int iter = 0; iter < 3; ++iter) { - sum_w = 0.0; sum_wI = 0.0; + sum_w = 0.0; sum_wI = 0.0; sum_wb = 0.0; sum_cwb = 0.0; for (size_t m = i; m < kk; ++m) { const auto &r2 = partials[ev[m]]; const double corr = r2.corr; const double I_corr = pooled_I(r2) * corr; const double sigma_corr = static_cast(r2.sigma) * 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. + // The non-signal variance as the integrator measured it, which is where it belongs: + // the box sum knows it directly and the profile fit knows the variance of its own + // fit, and neither is recoverable here once a sigma has been 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 a_var = bkg_var > 0.0 ? bkg_var : sigma_corr * sigma_corr; + double var = a_var + corr * std::max(0.0, F); const double w = 1.0 / var; sum_w += w; sum_wI += w * I_corr; + sum_wb += 1.0 / a_var; + sum_cwb += corr / (a_var * a_var); } F = sum_wI / sum_w; } + const double var_bkg_full = 1.0 / sum_wb; const int n_frames_event = static_cast(kk - i); i = kk; @@ -1608,6 +1613,8 @@ void RotationScaleMerge::Combine() { full.h = rawrun_h[r]; full.k = rawrun_k[r]; full.l = rawrun_l[r]; full.I = static_cast(F); full.sigma = static_cast(sigma_full); + full.var_bkg = static_cast(var_bkg_full); + full.var_per_I = static_cast(var_bkg_full * var_bkg_full * sum_cwb); full.d = d; full.rlp = 1.0f; full.partiality = 1.0f; @@ -2028,10 +2035,20 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool }; estimate_asymptote(); - auto corrected_sigma = [&](float I_corr, float sigma_corr, int g) -> float { + auto corrected_sigma = [&](const Obs &o, float I_corr, float sigma_corr) -> float { if (!error_model_active) return sigma_corr; - const double I_for_b = std::isfinite(em_mean[g]) ? em_mean[g] : I_corr; - const double v = error_model_a * static_cast(sigma_corr) * sigma_corr + const double I_for_b = std::isfinite(em_mean[o.group]) ? em_mean[o.group] : I_corr; + // A full's own sigma carries its own Poisson fluctuation: a full that came out low has a + // smaller sigma, so 1/sigma^2 weights it up and the merged mean drifts below . Rebuild the + // variance at the reflection's EXPECTED intensity - var_bkg + var_per_I*, the linear model + // the combine measured - so the weight no longer knows this full's own fluctuation. Mirrors + // MergeOnTheFly::CorrectedSigma on the stills path. + double a_var = static_cast(sigma_corr) * sigma_corr; + const double base = static_cast(o.corr) * o.corr * o.var_bkg + + static_cast(o.corr) * o.var_per_I * std::max(0.0, I_for_b); + if (base > 0.0) + a_var = base; + const double v = error_model_a * a_var + (error_model_b * I_for_b) * (error_model_b * I_for_b); return v > 0.0 ? static_cast(std::sqrt(v)) : sigma_corr; }; @@ -2077,7 +2094,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (!usable_merge(o)) continue; const float I_corr = o.I * o.corr; float sigma_corr = o.sigma * o.corr; - sigma_corr = corrected_sigma(I_corr, sigma_corr, o.group); + sigma_corr = corrected_sigma(o, I_corr, sigma_corr); if (reject_outliers && error_model_active && std::isfinite(reject_median[o.group]) && std::fabs(I_corr - reject_median[o.group]) > reject_nsigma * sigma_corr) { ++reject_count; @@ -2316,7 +2333,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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; - const float sigma_corr = corrected_sigma(I_corr, o.sigma * o.corr, o.group); + const float sigma_corr = corrected_sigma(o, I_corr, o.sigma * o.corr); if (!(sigma_corr > 0.0f) || !std::isfinite(sigma_corr)) continue; const double w = 1.0 / (static_cast(sigma_corr) * sigma_corr); AnomAcc &a = anom[HKLKey{ak.h, ak.k, ak.l, true}.pack()]; @@ -2631,15 +2648,18 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { fulls.assign(nf, Obs{}); std::vector fh(nf), fk(nf), fl(nf), fframe(nf), fgroup(nf); std::vector fI(nf), fsig(nf), fd(nf), fimg(nf), fcorr(nf, 1.0f), fpx(nf), fpy(nf); + std::vector fvb(nf), fvi(nf); std::vector fon(nf); gpu_->GetFulls(fh.data(), fk.data(), fl.data(), fI.data(), fsig.data(), fd.data(), fimg.data(), fframe.data(), fon.data(), fgroup.data()); gpu_->GetFullsPxPy(fpx.data(), fpy.data()); + gpu_->GetFullsVariance(fvb.data(), fvi.data()); if (scaled_fulls_on_gpu) gpu_->GetFullsCorr(fcorr.data()); for (int i = 0; i < nf; ++i) { Obs &o = fulls[i]; o.h = fh[i]; o.k = fk[i]; o.l = fl[i]; o.I = fI[i]; o.sigma = fsig[i]; o.d = fd[i]; + o.var_bkg = fvb[i]; o.var_per_I = fvi[i]; o.rlp = 1.0f; o.partiality = 1.0f; o.corr = fcorr[i]; o.image_number = fimg[i]; o.frame = fframe[i]; o.px = fpx[i]; o.py = fpy[i]; o.on_ice = fon[i]; o.group = fgroup[i]; diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 86f64461..242b6fa6 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -88,6 +88,9 @@ private: struct Obs { int32_t h, k, l; float I, sigma, d, rlp, partiality, zeta, delta_phi, bkg, var_bkg; + // Fulls only, written by the combine: the full's variance as a function of intensity, + // var(I) = var_bkg + var_per_I * I. The merge rebuilds it at the reflection's mean. + float var_per_I = 0.0f; 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 7f46bbfa..4c3166c8 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -266,7 +266,7 @@ namespace { int32_t *rr_nevents, *rr_nusable; // count pass outputs const int32_t *rr_offset; // emit pass: per-run base offset into the fulls arrays int32_t *f_h, *f_k, *f_l, *f_frame, *f_group; - float *f_I, *f_sigma, *f_d, *f_img, *f_px, *f_py; + float *f_I, *f_sigma, *f_d, *f_img, *f_px, *f_py, *f_var_bkg, *f_var_per_I; uint8_t *f_on_ice; }; @@ -346,27 +346,30 @@ namespace { } double F = sum_wI / sum_w; - // Pass C: 3 de-biased Poisson reweights (variance = bkg part + corr*max(0,F)). + // Pass C: 3 de-biased Poisson reweights (variance = bkg part + corr*max(0,F)), plus the + // full's var(I) = var_bkg + var_per_I*I for the merge (see the host combine). + double sum_wb = 0.0, sum_cwb = 0.0; for (int iter = 0; iter < 3; ++iter) { - sum_w = 0.0; sum_wI = 0.0; + sum_w = 0.0; sum_wI = 0.0; sum_wb = 0.0; sum_cwb = 0.0; for (int m = ev_start; m <= ev_end; ++m) { const int i = p.perm[m]; if (!CombineUsable(i, p.I, p.sigma, p.corr)) continue; const double corr = p.corr[i]; const double I_corr = pooled_I(i) * corr; 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. // 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 a_var = bkg_var > 0.0 ? bkg_var : sigma_corr * sigma_corr; + double var = a_var + corr * Dmax(0.0, F); const double w = 1.0 / var; sum_w += w; sum_wI += w * I_corr; + sum_wb += 1.0 / a_var; + sum_cwb += corr / (a_var * a_var); } F = sum_wI / sum_w; } + const double var_bkg_full = 1.0 / sum_wb; if (sum_w <= 0.0 || sum_partiality < p.min_partiality || sum_partiality < p.min_captured_fraction) @@ -384,6 +387,8 @@ namespace { p.f_h[o] = p.rr_h[r]; p.f_k[o] = p.rr_k[r]; p.f_l[o] = p.rr_l[r]; p.f_I[o] = float(F); p.f_sigma[o] = float(sigma_full); + p.f_var_bkg[o] = float(var_bkg_full); + p.f_var_per_I[o] = float(var_bkg_full * var_bkg_full * sum_cwb); p.f_d[o] = d; p.f_img[o] = peak_frame; p.f_px[o] = peak_px; p.f_py[o] = peak_py; @@ -429,7 +434,7 @@ namespace { int n_groups; double min_partiality, error_model_a, error_model_b, reject_nsigma; int for_search, error_model_active, reject_outliers; - const float *I, *sigma, *corr, *partiality, *d, *reject_median; + const float *I, *sigma, *corr, *partiality, *d, *reject_median, *var_bkg, *var_per_I; const int32_t *group, *frame; const uint8_t *on_ice, *frame_cell_ok; const int32_t *gperm, *gstart, *gcount; @@ -523,7 +528,13 @@ namespace { if (p.error_model_active) { const double I_for_b = isfinite(p.em_mean[g]) ? p.em_mean[g] : double(I_corr); const double bi = p.error_model_b * I_for_b; - const double v = p.error_model_a * double(sigma_corr) * sigma_corr + bi * bi; + // The full's variance rebuilt at the reflection's expected intensity, so the merge + // weight cannot know this full's own fluctuation (see the host corrected_sigma). + const double c = p.corr[i]; + double a_var = double(sigma_corr) * sigma_corr; + const double base = c * c * p.var_bkg[i] + c * p.var_per_I[i] * Dmax(0.0, I_for_b); + if (base > 0.0) a_var = base; + const double v = p.error_model_a * a_var + bi * bi; if (v > 0.0) sigma_corr = float(sqrt(v)); } if (p.reject_outliers && p.error_model_active && isfinite(rmed) @@ -627,7 +638,7 @@ struct RotationScaleMergeGPU::Impl { // combine: resident fulls SoA (rebuilt each Combine) int n_fulls = 0; CudaDevicePtr f_h, f_k, f_l, f_frame, f_group; - CudaDevicePtr f_I, f_sigma, f_d, f_img, f_px, f_py; + CudaDevicePtr f_I, f_sigma, f_d, f_img, f_px, f_py, f_var_bkg, f_var_per_I; CudaDevicePtr f_on_ice; // scale-fulls (Unity model, kept resident): all-ones partiality/rlp/zeta so the shared scaling kernels // yield coeff=mean, plus the working corr, the per-obs scale scratch, and the fulls frame/group CSRs @@ -794,6 +805,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); + p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.em_mean = d.m_em_mean.get(); p.a_swI = d.a_swI.get(); p.a_sw = d.a_sw.get(); p.a_swIh0 = d.a_swIh0.get(); p.a_swIh1 = d.a_swIh1.get(); @@ -944,6 +956,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti d.f_I = CudaDevicePtr(nf); d.f_sigma = CudaDevicePtr(nf); d.f_d = CudaDevicePtr(nf); d.f_img = CudaDevicePtr(nf); d.f_px = CudaDevicePtr(nf); d.f_py = CudaDevicePtr(nf); + d.f_var_bkg = CudaDevicePtr(nf); d.f_var_per_I = CudaDevicePtr(nf); d.f_on_ice = CudaDevicePtr(nf); d.f_corr = CudaDevicePtr(nf); d.f_partiality = CudaDevicePtr(nf); d.f_rlp = CudaDevicePtr(nf); d.f_zeta = CudaDevicePtr(nf); @@ -956,6 +969,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.f_frame = d.f_frame.get(); p.f_group = d.f_group.get(); p.f_I = d.f_I.get(); p.f_sigma = d.f_sigma.get(); p.f_d = d.f_d.get(); p.f_img = d.f_img.get(); p.f_px = d.f_px.get(); p.f_py = d.f_py.get(); + 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(); if (d.n_fulls > 0) { CombineKernel<<>>(p); @@ -1053,6 +1067,15 @@ void RotationScaleMergeGPU::GetFullsPxPy(float *px, float *py) const { CudaCheck(cudaMemcpy(py, d.f_py.get(), bytes, cudaMemcpyDeviceToHost), "download f_py"); } +void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I) const { + 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), + "download f_var_per_I"); +} + void RotationScaleMergeGPU::SetFullsCorr(const float *corr) { auto &d = *impl_; if (d.n_fulls == 0) return; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 9aea3739..337e737b 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -125,6 +125,9 @@ public: // surface. Length = n_fulls. void GetFullsPxPy(float *px, float *py) const; + // Download the fulls' variance model, var(I) = var_bkg + var_per_I * I. Length = n_fulls. + void GetFullsVariance(float *var_bkg, float *var_per_I) const; + // Re-upload the fulls' working corr (length n_fulls) after the host correction surfaces (decay / // absorption) mutate it, so the resident merge reads the corrected scale. void SetFullsCorr(const float *corr);