From d8029524e76f456c0be507fe2407c53034b915a1 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 11 Aug 2026 16:06:39 +0200 Subject: [PATCH] Scaling: never let an observation's own fluctuation set its weight A weighted mean is only unbiased while the weights are independent of the values being averaged. The IUCr's own nomenclature report (Schwarzenbach et al., Acta Cryst A45 (1989) 63-75) puts it directly: weights in averaging "should not be based on the counting statistics of the individual observations whose estimated variances are biased and result in larger weights for accidentally low intensities". Two places in the rotation pipeline were doing exactly that, and between them they drove whole resolution shells of merged intensity negative. 1. The profile fit computed its non-signal variance as var_bkg = max(0, 1/den - max(0, I) + bkg-estimate term) The point of a separate var_bkg is that it does NOT move with the reflection's own fluctuation, and 1/den - I is the quantity that does not: 1/den is the fit variance taken at the fitted intensity and grows with it roughly one for one. Clamping the subtrahend at zero left a down-fluctuated reflection's own deflated variance standing as its background variance. Measured over 6.9 M partials of one weak rotation dataset, var_bkg/bkg came out at 3.7-5.4 for observations with I < 0 against 11.4-13.7 for I > 0 - the down-fluctuated half of every reflection carried a variance ~2.7x too small and was weighted up by the same factor, first in the 3D combine and then again in the merge. Removing the clamp makes var_bkg flat in I (~13 x bkg across the whole range). 2. The merge then weighted each combined full by 1/sigma_full^2, and sigma_full is by construction a function of the full's own answer: the combine's variance carries a corr*max(0, F) signal term, so every full with F <= 0 got the smallest variance the model allows while the strongest quartile got 2.26x more. The merge now rebuilds that variance at the reflection's mean instead, from a linear model var(I) = var_bkg + var_per_I * I that the combine measures and stores on the full. This mirrors MergeOnTheFly::CorrectedSigma, whose comment already claimed to mirror the rotation combine. Verified against an estimator that cannot see the fluctuation - summing the partials and dividing by the summed partiality, the classical construction every other program uses (Greenhough & Suddath, J. Appl. Cryst. 19 (1986) 400-409, via Leslie, Acta Cryst D55 (1999) 1696-1702: profile fitting biases the individual partials but not their sum). Reproducing the merge on dumped observations, the shipped weighting sat ~1.9 sigma below that reference in the noise shells; the two changes recover most of it, and every intensity-independent weighting scheme agrees with the reference once (1) is in. Four-crystal probe, XDS resolution limits, branch fingerprint identical on all four (so none of these is a two-pass branch flip): weak cubic case last shell -1.6 -> +0.2 (XDS +0.10), last shell R_meas 478% -> 250% (XDS 246%), overall 6.1 -> 7.5 (XDS 7.18), R_meas 18.3% -> 18.1%, CC1/2_hi 38.2% -> 43.7% tetragonal case outer shells -0.4/-0.8/-0.9/-1.0 -> +1.8/+1.2/ +0.9/+0.4, R_meas 184%/595%/7614%/nan -> 95%/119%/135%/232% (the nan was the shell mean crossing zero), R_meas 33.3% -> 32.9%, CC1/2_hi 38.3% -> 56.5% trigonal case R_meas 13.0% -> 12.5%, CC1/2_hi 14.4% -> 16.5% strong control unchanged to every printed digit but ISa Cost: ISa falls (17.2 -> 14.0 and 16.7 -> 14.9 on the two mid-strength cases, 28.3 -> 27.8 on the control). Strong reflections are untouched by (1) - their partials are all positive, so var_bkg is bit-identical - but the joint a/b fit redistributes: honest weak sigmas lower a, and b rises to keep the strong bins fitted. The median reduced chi^2 improves (1.25 -> 1.14, 1.35 -> 1.28) so the new split describes the scatter better, but ISa is the one headline metric that moves the wrong way and it should be watched over the full battery. The integrator change is shared, so the stills merge sees it too; there it feeds GetExpectedVarianceMerge, which had been handed the same contaminated var_bkg. That path is untested here. Co-Authored-By: Claude Opus 5 (1M context) --- docs/CHANGELOG.md | 2 + .../BraggIntegrationEngineCPU.cpp | 8 +++- .../BraggIntegrationEngineGPU.cu | 3 +- .../scale_merge/RotationScaleMerge.cpp | 44 ++++++++++++++----- .../scale_merge/RotationScaleMerge.h | 3 ++ .../scale_merge/RotationScaleMergeGPU.cu | 43 +++++++++++++----- .../scale_merge/RotationScaleMergeGPU.h | 3 ++ 7 files changed, 82 insertions(+), 24 deletions(-) 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);