mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 68 Wiki Activity

Scaling: fit a correction surface with the observation as the response, not the regressor

Browse Source 
ApplyCellSurface fits one multiplicative factor per cell by least squares with the
OBSERVATION on the regressor side - A = sum w Is Iref / sum w Is^2, the slope that
carries Is onto Iref. A least-squares slope is attenuated by the noise in its own
regressor, here by 1/(1 + (sigma/I)^2), and an observation carries all of a
reflection's noise where the reference carries about 1/n of it. So every cell is
pulled towards zero by an amount set by its own signal-to-noise - and on a detector
that is a function of radius, which is to say of resolution. The gauge fix then
spreads the ramp over the whole surface and the alternating rounds compound it.

Nothing downstream catches it. The cross-validation splits by frame parity, and a
bias that depends only on a cell's signal-to-noise is identical in both halves. And
the held-out score is sum|Is - Iref| / sum Iref over the whole resolution range,
which any resolution-dependent scale lowers without tightening a single reflection:
applying a scale that is purely a function of d to a merged 360 deg sweep leaves
every resolution shell's R_meas unchanged to 0.05 pp and takes the run's overall
R_meas from 42.4% to 27.5%. The surface finds that manipulation because its own bias
points exactly along it, and reports it as a 40% held-out gain.

The cost is large wherever a sweep was taken past its signal. On one such run the
fitted detector-plane "flat field" ran from 0.25 to 4.0 with 37% of its cells pinned
at the low clamp - an 11x centre-to-edge ramp - and on the same combined fulls it
moved the merge: low-resolution R_meas 8.8% -> 15.5%, <I/sigma> 33.2 -> 13.4, CC1/2
99.65 -> 98.87, error model b 5.0e-03 -> 4.3e-02, ISa 11.0 -> 3.5. The program's own
--no-scaling-corrections run agrees on the same fulls and the same space group
(b 5.4e-03 ISa 10.6 against b 4.1e-02 ISa 3.6). Regressing the other way round
restores all of it - 8.6%, 33.6, 99.65, 4.9e-03, ISa 11.1 - and keeps the surface's
real gain in the middle shells, where CC1/2 rises by 1-2 points.

Where the data are well measured everywhere the two estimators are indistinguishable:
on a control crystal the two surfaces agree to 0.02 pp in every shell and 0.02 in
ISa, which is what a flat field should look like. Recovering a synthetic +-20% ripple
imposed on the same fulls: 0.18 rms in log against 0.91 for the shipped form on the
weak-outer-shell crystal, 0.059 against 0.082 on the well-measured one, with the
spurious correlation between the fitted factor and detector radius falling from -0.30
to -0.01.

The same estimator serves the absorption surface and any resolution-indexed surface
fitted through this function, where the bias lands directly on the resolution axis:
on a crystal the program itself reports as having no radiation damage (total dB
0.00 A^2), a batch x resolution-shell surface fitted the old way manufactures a
monotone 12% falloff from low to high resolution out of nothing, and the new way
gives 1%.

The per-frame scale (FitPerFrameG) already regresses this way round.

Full 38-crystal rotation battery against the same binary without it. Read it with
the next commit, which supplies the correction this one stops faking; alone it
is a partial state:

  observations  better 34 / worse  4,  +22 940
  R_meas_lo     better  4 / worse  8,     -5.9
  ISa           better 20 / worse 14,   +32.66
  CC1/2         better  2 / worse 10,    -11.4
  R_meas        better  1 / worse 30,   +114.3
  space groups  unchanged at 35/38

Overall R_meas RISES, and that is the artefact leaving rather than arriving: it
is a ratio of sums across every shell, so a resolution-dependent scale lowers it
without one reflection getting tighter - measured, a scale that is purely f(d)
leaves every shell's R_meas unchanged to 0.05 pp while moving the run's overall
value from 42.4 to 27.5 per cent. That is precisely the shape of the bias, which
is why the surface's own held-out score read it as a 40 per cent gain.

CC1/2 falling on ten crystals is not covered by that argument and is the reason
this commit is not defensible on its own: the ramp was partly standing in for a
real time-dependent absorption that nothing else modelled. With that correction
supplied by the following commit the same battery gives CC1/2 -1.4 and
R_meas_lo -26.0.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-08-12 05:15:49 +02:00
co-authored by Claude Opus 5
parent d74ff7fc29
commit b8e7e9c8bf
2 changed files with 17 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+16 -6
View File
@@ -1567,7 +1567,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
// Fit the per-cell factor over the subset {frame&1 == parity} (parity < 0 = all fulls), n_iter
// alternating rounds against that subset's own reference (Tikhonov pull to 1, gauge-fixed to a
// den-weighted geometric mean of 1 so it never drifts the overall scale). Returns the per-cell factor.
// weighted geometric mean of 1 so it never drifts the overall scale). Returns the per-cell factor.
auto fit_surface = [&](int parity) -> std::vector<double> {
const std::vector<int32_t> &sel = subset(parity);
std::vector<double> A(ncell, 1.0);
@@ -1579,7 +1579,17 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
const double a = A[cell[i]], sc = static_cast<double>(o.sigma) * o.corr * a, w = 1.0 / (sc * sc);
sw[o.group] += w; swI[o.group] += w * static_cast<double>(o.I) * o.corr * a;
}
std::vector<double> num(ncell, 0.0), den(ncell, 0.0);
// The cell's factor is the scale that carries the observation onto the reference, and it is
// fitted by regressing the OBSERVATION on the reference: A = sum w Iref^2 / sum w Is Iref,
// the inverse of the slope of Is against Iref. Not the other way round - a least-squares
// slope is attenuated by the noise in its own regressor, by 1/(1 + (sigma/I)^2), and an
// observation carries all of a reflection's noise where the reference carries about 1/n of
// it. That attenuation is a property of the cell's signal-to-noise, hence of where the cell
// sits on the detector, so with Is as the regressor the bias is a smooth radial ramp - which
// the gauge fix spreads over the whole surface and the n_iter rounds compound, and which no
// even/odd cross-validation can see because both halves carry it equally. The per-frame
// scale (FitPerFrameG) already regresses this way round.
std::vector<double> cross(ncell, 0.0), ref2(ncell, 0.0);
for (const int32_t i : sel) {
const Obs &o = fulls[i];
if (sw[o.group] <= 0.0) continue;
@@ -1587,15 +1597,15 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
const double Is = static_cast<double>(o.I) * o.corr * a, sc = static_cast<double>(o.sigma) * o.corr * a;
if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0) || !(sc > 0.0)) continue;
const double w = 1.0 / (sc * sc);
num[cell[i]] += w * Is * Iref; den[cell[i]] += w * Is * Is;
cross[cell[i]] += w * Is * Iref; ref2[cell[i]] += w * Iref * Iref;
}
std::vector<double> dsorted = den;
std::vector<double> 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<double> upd(ncell, 1.0);
for (int c = 0; c < ncell; ++c) upd[c] = (num[c] + lambda) / (den[c] + lambda);
for (int c = 0; c < ncell; ++c) upd[c] = (ref2[c] + lambda) / (cross[c] + lambda);
double logsum = 0.0, wsum = 0.0;
for (int c = 0; c < ncell; ++c) if (den[c] > 0.0) { logsum += den[c] * std::log(upd[c]); wsum += den[c]; }
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;
for (int c = 0; c < ncell; ++c) A[c] = std::clamp(A[c] * upd[c] / gm, 0.25, 4.0);
}