From 20349fb2d97ff4764a584f42ba5ef3b46ae562b2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 21 Sep 2026 17:11:44 +0200 Subject: [PATCH] Scaling: a correction surface is fitted under a prior and gated on held-out half-set CC Three changes to ApplyCellSurface, each measured offline on the surfaces' own input terms from twelve open rotation sets (36 surfaces), against per-shell CC of the merge to the deposited model's Fc. The gate. The held-out score was sum|Is - Iref| / sum Iref pooled over every shell the fit reaches. The surfaces are fitted before the resolution cutoff, so on a weak set the score was carried by shells with no signal (one set: 8 of 20 shells at < 0.1, a quarter of the numerator), where any spread in A raises |Is - Iref| because Is is A times noise; and it weighted each residual equally where the merge weights it by 1/(sigma corr A)^2. It is replaced by the quantity the surface is for: fit on even frames and merge the odd ones with it, fit on odd and merge the even ones, and correlate the two half-set means within each resolution shell, against the same with no surface. A surface fitted to noise lowers that correlation, so the gate needs no margin. On surfaces moving the external CC by more than 0.005 either way, the old gate agreed with the model on 6 of 12 and the new one on 13 of 14; and no threshold on the old score separated them - below its 2% sat surfaces worth +0.017 and +0.009 and surfaces worth -0.026 and -0.013. CV_MIN_RELATIVE_GAIN still gates the decay and relative-B fits, which this does not touch. The regulariser. The per-round update A *= (ref2 + lambda) / (cross + lambda) was described as a Tikhonov pull to 1, but lambda cancels at the fixed point ref2 = cross: it only slowed a weak cell to cross / (cross + lambda) of its step, so the 30-round cap regularised by early stopping and the weak cells were still moving when it hit - most of what "NOT settled" reported - and a cell with no information was fitted to its noise and pinned on the [0.25, 4] clamp. The update is now a damped Gauss-Newton step on log A with a N(0, 0.1^2) prior: ref2 is the information of the cell's log scale, so the fixed point ref2 - cross = P log A moves a cell off 1 in proportion to what its observations say. Offline, applied surfaces that settle within the cap: most of them, against a minority; cells on the clamp: 0-4 per surface against up to 318. On the weakest set the time x detector surface goes from 181 clamped cells and -0.031 external CC to none and +0.001. The gauge. The per-shell projection is exact only where a cell sits at one resolution - a detector patch. A goniometer-frame cell collects directions from the whole sweep, so there it removed real absorption, which varies with scattering angle: +0.016 -> +0.022 and +0.017 -> +0.020 on the two sets where that surface is worth anything. That surface keeps only the global level fix. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- .../scale_merge/RotationScaleMerge.cpp | 157 +++++++++++------- .../scale_merge/RotationScaleMerge.h | 8 +- 2 files changed, 103 insertions(+), 62 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 67a8582b8..4e7be9f8c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -2973,7 +2973,8 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) { cell[i] = ix * NB + iy; } }); - ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Absorption (goniometer-frame 8x8)"); + ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Absorption (goniometer-frame 8x8)", + /*resolution_gauge=*/false); } void RotationScaleMerge::RefineAbsorptionTime(int n_iter, int n_groups) { @@ -3053,7 +3054,8 @@ void RotationScaleMerge::RefineAbsorptionTime(int n_iter, int n_groups) { cell[i] = (it * ND + ix) * ND + iy; } }); - ApplyCellSurface(cell, NT * ND * ND, n_iter, n_groups, "Absorption (time x detector 12x10x10)"); + ApplyCellSurface(cell, NT * ND * ND, n_iter, n_groups, "Absorption (time x detector 12x10x10)", + /*resolution_gauge=*/true); } void RotationScaleMerge::RefineModulation(int n_iter, int n_groups) { @@ -3091,11 +3093,12 @@ void RotationScaleMerge::RefineModulation(int n_iter, int n_groups) { cell[i] = ix * NB + iy; } }); - ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Modulation (detector-frame 16x16)"); + ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Modulation (detector-frame 16x16)", + /*resolution_gauge=*/true); } void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int ncell, int n_iter, - int n_groups, const char *name) { + int n_groups, const char *name, bool resolution_gauge) { auto usable = [&](const Obs &o) { return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && o.partiality >= min_partiality; }; @@ -3188,6 +3191,14 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int // same simulation that recovers a 12% cos(2 phi) flat field to CC 0.986 at full amplitude // while the radial ramp goes to 1.00. The projection removes the overall level too, so it // subsumes the global gauge below. + // It is only a null direction where a cell sits at one resolution - a patch of the detector, which + // is what the modulation and the time x detector cells are. A goniometer-frame cell collects + // diffracted directions from the whole sweep, so its observations span most of the resolution + // range; the shell means are then a mix over cells rather than a direction the data cannot see, + // and projecting them out removes real absorption, which does vary with scattering angle. Measured + // on four open sets (external per-shell CC against the deposited model, the surface applied alone): + // +0.016 -> +0.022 and +0.017 -> +0.020 on the two where the surface is worth anything, no change on + // the other two. That surface keeps the global level fix only (resolution_gauge false). // Equal-occupancy shells in s^2 = 1/d^2. The count is second order (8, 20 and 40 shells all leave // the residual ramp inside 1.5%); 20 is fine enough to follow a curved drift and coarse enough // that each shell still holds tens of thousands of observations. @@ -3198,6 +3209,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int // the alternating fit below converges against it. std::vector shw(static_cast(ncell) * nshell, 0.0); std::vector shw_shell(nshell, 0.0), shw_cell(nshell > 0 ? ncell : 0, 0.0); + std::vector group_shell(n_groups, 0); // the shell each ASU group sits in, for the gate if (nshell > 0) { std::vector s2; s2.reserve(term.size()); @@ -3215,6 +3227,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int const Term &t = term[k]; const float v = t.d > 0.0f ? 1.0f / (t.d * t.d) : 0.0f; const int s = static_cast(std::upper_bound(edge.begin(), edge.end(), v) - edge.begin()); + group_shell[t.group] = s; const double sc = static_cast(t.sigma) * t.corr; if (!(sc > 0.0)) continue; shw[static_cast(t.cell) * nshell + s] += 1.0 / (sc * sc); @@ -3269,8 +3282,8 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int }; // Fit the per-cell factor over the subset {frame&1 == parity} (parity < 0 = all fulls), alternating - // rounds against that subset's own reference (Tikhonov pull to 1, then gauge-fixed within each - // resolution shell, which also holds the overall scale) until the factors settle - the same + // rounds against that subset's own reference (a prior pull to 1, then gauge-fixed within each + // resolution shell where resolution_gauge, which also holds the overall scale) until the factors settle - the same // tolerance as the per-frame scales, on the rms |log| of the round's update - or n_iter rounds // have been spent. Returns the per-cell factor. int rounds_used = 0; @@ -3335,21 +3348,36 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int std::vector dsorted = cross; std::nth_element(dsorted.begin(), dsorted.begin() + dsorted.size() / 2, dsorted.end()); const double lambda = 0.1 * std::max(1e-30, dsorted[dsorted.size() / 2]); - std::vector upd(ncell, 1.0); - // A cell whose observations carry no net signal can sum to `cross` <= 0 now that the - // negatives are in, and the Tikhonov pull does not always lift it clear of zero. Such a - // cell has nothing to say about its own factor, so it keeps the one it has rather than - // dividing by a number that can be near zero or negative - which would reach the log below - // as a NaN or as a factor pinned on the clamp and, through the shell means, poison every - // other cell. It is the same `cross > 0` the gauge and the step size already test. - for (int c = 0; c < ncell; ++c) - if (cross[c] > 0.0) upd[c] = (ref2[c] + lambda) / (cross[c] + lambda); - double logsum = 0.0, wsum = 0.0; - for (int c = 0; c < ncell; ++c) if (cross[c] > 0.0) { logsum += cross[c] * std::log(upd[c]); wsum += cross[c]; } - const double gm = wsum > 0.0 ? std::exp(logsum / wsum) : 1.0; + // One damped Gauss-Newton step on log A per round, with a N(0, SURFACE_PRIOR_SIGMA^2) prior on + // each cell's log factor. ref2 is the Fisher information of the cell's log scale (the inverse + // variance of the slope of Is against Iref) and ref2 - cross is its first-order residual, of + // either sign - so a cell whose negatives outweigh its positives needs no special case. + // The prior is what regularises. `lambda` used to be called a Tikhonov pull, but in the old + // update (ref2 + lambda) / (cross + lambda) it cancels at the fixed point, ref2 = cross: it only + // slowed a weak cell to cross / (cross + lambda) of its step, so the 30-round cap did the + // regularising as early stopping, and a weak cell was still walking when the cap hit it - + // most of what "NOT settled" reported. Here lambda is kept for that damping alone (the + // undamped step oscillates on some surfaces), and the prior sets the fixed point: + // ref2 - cross = P log A, so a cell's factor moves off 1 in proportion to the information it + // carries against the fixed prior precision P - a cell whose observations say little stays + // near 1 instead of being fitted to its noise and pinned on the clamp. + // SURFACE_PRIOR_SIGMA = 0.1 against 0.2 on six open sets, external per-shell CC against the + // deposited model: 0.1 is as good or better on 13 of 18 surfaces, and settles in fewer rounds. + constexpr double SURFACE_PRIOR_SIGMA = 0.1; + const double P = 1.0 / (SURFACE_PRIOR_SIGMA * SURFACE_PRIOR_SIGMA); std::vector lnext(ncell); - for (int c = 0; c < ncell; ++c) lnext[c] = std::log(A[c] * upd[c] / gm); - if (nshell > 0) { + for (int c = 0; c < ncell; ++c) { + lnext[c] = std::log(A[c]); + if (ref2[c] > 0.0) + lnext[c] += (ref2[c] - cross[c] - P * lnext[c]) / (ref2[c] + P + lambda); + } + if (!(resolution_gauge && nshell > 0)) { + // The overall level is a null direction of every surface: fix it at 1. + double acc = 0.0, wsum = 0.0; + for (int c = 0; c < ncell; ++c) if (ref2[c] > 0.0) { acc += ref2[c] * lnext[c]; wsum += ref2[c]; } + if (wsum > 0.0) + for (int c = 0; c < ncell; ++c) lnext[c] -= acc / wsum; + } else { // Project out the component that is a function of resolution: the weighted mean of // log A within each shell, carried back to a cell in the proportion in which that // cell's observations sit in each shell. @@ -3373,7 +3401,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int n_clamped = 0; for (int c = 0; c < ncell; ++c) { const double next = std::clamp(std::exp(lnext[c]), 0.25, 4.0); - if (cross[c] > 0.0) { + if (ref2[c] > 0.0) { step2 += std::pow(std::log(next / A[c]), 2); ++n_step; if (next <= 0.2500001 || next >= 3.9999999) ++n_clamped; @@ -3385,50 +3413,61 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int break; } } - if (nshell > 0) + if (resolution_gauge && nshell > 0) gauge_ramp = std::exp(*std::max_element(removed.begin(), removed.end()) - *std::min_element(removed.begin(), removed.end())); return A; }; - // Mean studentized squared deviation over the subset {frame&1 == parity} with surface A applied, - // scored against that subset's own reference (no leakage). Lower = tighter equivalents. - auto score = [&](int parity, const std::vector &A) -> double { - const std::vector &sel = subset(parity); + // Cross-validate on the half-set agreement it is meant to improve. Fit the surface on the even + // frames and merge the ODD frames with it; fit it on the odd frames and merge the EVEN frames with + // it; then correlate the two half-set means within each resolution shell, and compare with the same + // two halves merged with no surface. Each half is corrected by a surface it did not help to fit, so + // a surface fitted to noise adds independent noise to both and lowers the correlation, and only a + // correction the two halves share raises it: no margin is needed, the null already sits below zero. + // Within a shell a correlation is blind to a scale that is a function of resolution, which is the + // one component this surface cannot determine. + // It replaces a held-out sum|Is - Iref| / sum Iref over the whole range, which two things made + // disagree with the external referent. Pooled over every shell the fit reaches, it was carried by + // shells past the resolution the merge keeps: on one set 8 of 20 shells held < 0.1 and a + // quarter of the numerator, where any spread in A raises |Is - Iref| (Is is A times noise) whether + // the surface is right or not. And it counted each residual at full weight where the merge weights it + // by 1/(sigma corr A)^2. Measured on twelve open sets, each surface alone against per-shell CC to the + // deposited model: of the surfaces that move that CC by more than 0.005 either way, the old gate + // took the model's side on 6 of 12 and this one on 13 of 14. The old gate's 2% margin was not what + // stood between the two: below it sat surfaces worth +0.017 and +0.009 and surfaces worth -0.026 + // and -0.013 alike, so no threshold on that score separates them. + const int nsh_cc = nshell > 0 ? nshell : 1; + auto half_means = [&](int parity, const std::vector &A) { reference(parity, A); - // Rmeas-like (sigma-INDEPENDENT) agreement of the held-out equivalents: sum|Is - Iref| / sum|Iref|. - // Scoring on the studentized deviation instead lets a surface "improve" the held-out chi^2 by - // reshaping sigma (via corr) without tightening the actual intensities - which on mis-indexed / bad - // data passes cross-validation yet worsens Rmeas. A fractional metric cannot be gamed that way. - std::vector tnum(nt, 0.0), tden(nt, 0.0); - ParallelChunks(nt, nthreads, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - double num = 0.0, den = 0.0; - for (size_t k = part(sel.size(), t); k < part(sel.size(), t + 1); ++k) { - const Term &o = term[sel[k]]; - if (sw[o.group] <= 0.0) continue; - const double a = A[o.cell], Is = static_cast(o.I) * o.corr * a; - const double Iref = swI[o.group] / sw[o.group]; - if (!std::isfinite(Iref) || Iref <= 0.0) continue; - num += std::abs(Is - Iref); den += Iref; - } - tnum[t] = num; tden[t] = den; - } - }); - double num = 0.0, den = 0.0; - for (int t = 0; t < nt; ++t) { num += tnum[t]; den += tden[t]; } - return den > 0.0 ? num / den : 0.0; + std::vector m(n_groups, std::numeric_limits::quiet_NaN()); + for (int g = 0; g < n_groups; ++g) + if (sw[g] > 0.0) m[g] = swI[g] / sw[g]; + return m; + }; + auto shell_cc = [&](const std::vector &x, const std::vector &y, int s) { + double n = 0, sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; + for (int g = 0; g < n_groups; ++g) { + if (group_shell[g] != s || !std::isfinite(x[g]) || !std::isfinite(y[g])) continue; + n += 1; sx += x[g]; sy += y[g]; sxx += x[g] * x[g]; syy += y[g] * y[g]; sxy += x[g] * y[g]; + } + const double vx = sxx - sx * sx / n, vy = syy - sy * sy / n; + return (n >= 50 && vx > 0.0 && vy > 0.0) ? (sxy - sx * sy / n) / std::sqrt(vx * vy) + : std::numeric_limits::quiet_NaN(); }; - - // Cross-validate: fit the surface on even frames and score the held-out odd equivalents (and vice - // versa). A real surface generalizes; an over-fit one (few obs / cell) does not. Apply only if the - // held-out equivalent chi^2 improves by a clear margin - so the correction can never worsen the data. const std::vector ident(ncell, 1.0); const std::vector A_even = fit_surface(0), A_odd = fit_surface(1); - const double base = score(1, ident) + score(0, ident); - const double gain = base - (score(1, A_even) + score(0, A_odd)); - if (!(gain > CV_MIN_RELATIVE_GAIN * base)) { - logger.Info("{} correction: not cross-validated (held-out gain {:.1f}%, skipped)", - name, 100.0 * gain / std::max(base, 1e-30)); + const std::vector odd0 = half_means(1, ident), even0 = half_means(0, ident); + const std::vector odd1 = half_means(1, A_even), even1 = half_means(0, A_odd); + double gain = 0.0; + int n_cc = 0; + for (int sh = 0; sh < nsh_cc; ++sh) { + const double c0 = shell_cc(odd0, even0, sh), c1 = shell_cc(odd1, even1, sh); + if (std::isfinite(c0) && std::isfinite(c1)) { gain += c1 - c0; ++n_cc; } + } + if (n_cc > 0) gain /= n_cc; + if (!(gain > 0.0)) { + logger.Info("{} correction: not cross-validated (held-out half-set CC {:+.4f} per shell, skipped)", + name, gain); return; } const std::vector A = fit_surface(-1); @@ -3437,9 +3476,9 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int if (cell[i] >= 0) fulls[i].corr = static_cast(fulls[i].corr * A[cell[i]]); }); - logger.Info("{} correction: cross-validated, held-out gain {:.1f}%; {} in {} round(s); " + logger.Info("{} correction: cross-validated, held-out half-set CC {:+.4f} per shell; {} in {} round(s); " "resolution gauge removed a {:.2f}x ramp; {} of {} cells on the [0.25, 4] clamp", - name, 100.0 * gain / std::max(base, 1e-30), settled ? "settled" : "NOT settled", + name, gain, settled ? "settled" : "NOT settled", rounds_used, gauge_ramp, n_clamped, ncell); } diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 68fa1bf83..03a4e75fe 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -582,10 +582,12 @@ private: // in the detector frame (not tied to the rotation) the same correction concept applies to stills. void RefineModulation(int n_iter, int n_groups); // Shared engine for the correction surfaces: given a per-full cell assignment (cell[i] in [0,ncell), or - // <0 to skip), fit a Tikhonov-regularised multiplicative factor per cell against the merged reference, - // cross-validate on even/odd frames, and fold it into corr only if the held-out equivalents improve. + // <0 to skip), fit a multiplicative factor per cell against the merged reference under a prior pull to 1, + // cross-validate on even/odd frames, and fold it into corr only if the held-out half-set CC improves. + // resolution_gauge: the cells are patches of the detector, each at one resolution, so the component of + // the surface that is a function of resolution is a null direction and is projected out per shell. void ApplyCellSurface(const std::vector &cell, int ncell, int n_iter, int n_groups, - const char *name); + const char *name, bool resolution_gauge); // Sort `fulls` by peak frame and (re)build fulls_frame_start/count (the per-frame CSR the scale-fulls // step slices). Shared by the CPU Combine tail and the GPU combine path.