diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b1fcd14d..50c4f58e 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.161 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. +* Scaling: the rotation correction surfaces (decay, absorption, detector-plane modulation) no longer pull each cell's factor towards zero in proportion to that cell's own noise, which had put a spurious centre-to-edge ramp on runs whose outer resolution shells carry no signal. * rugnux: A goniometer that turned further than the angles stored in the file — which are the commanded ones — is now measured and **corrected in the second rotation pass**, or set by hand with `--rotation-scale `; the fitted correction is applied only when it exceeds 0.5 %, moves each end of the sweep by at least 0.5°, and comes out the same with any fifth of the sweep left out. * rugnux: Every `mx` and `scale` run now writes `_report.txt`, a results report modelled on XDS's `CORRECT.LP` — `KEY= value` lines, fixed-width tables and `WARNING:` sentences covering indexing, geometry post-refinement, the space-group decision, merging, twinning, radiation damage and the stretches of the sweep over which the crystal delivered much less than the rest of the run. * 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. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index fe58bd68..08ce4a96 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1567,7 +1567,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &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 { const std::vector &sel = subset(parity); std::vector A(ncell, 1.0); @@ -1579,7 +1579,17 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int const double a = A[cell[i]], sc = static_cast(o.sigma) * o.corr * a, w = 1.0 / (sc * sc); sw[o.group] += w; swI[o.group] += w * static_cast(o.I) * o.corr * a; } - std::vector 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 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 &cell, int const double Is = static_cast(o.I) * o.corr * a, sc = static_cast(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 dsorted = den; + 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); - 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); }