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.