diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 10b0d0a9..2c2a8b2e 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -22,6 +22,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Packaging: Build the self-contained Linux viewer against a static libdbus with glib disabled, so the `.tgz` no longer drags a tail of system `.so` dependencies; add parallel image-build and in-container viewer-verification scripts. * rugnux: Write anomalous data as a standard CCP4 anomalous MTZ - one row per reflection with `IMEAN`, `I(+)`/`I(-)` (and the matching `F`/`F(+)`/`F(-)` amplitudes) - instead of two `IMEAN` rows per Bijvoet pair, so `aimless`/`mtz2sca`/ANODE read it directly. The non-anomalous MTZ output is unchanged. * rugnux: Always write the merged reflections as both `.mtz` and `.cif` (each has its uses downstream - MTZ for the CCP4/phenix tools, mmCIF for deposition). The `--scaling-output` format selector and the plain-text `.hkl` output are removed. +* rugnux: Add a detector-plane **modulation** (flat-field) correction surface to rotation scaling - a smooth multiplicative factor over the detector position where a reflection lands, fitted against the merged reference (symmetry-equivalents land at different positions as the crystal rotates) and cross-validated. It corrects detector-response/geometric systematics that inflate R-meas; on the rotation battery it lowers R-meas by several to tens of percent (e.g. lysoC 23->16%, lyso_2 47->28%, EcwtAL500 53->28%) with CC1/2 held or improved and the anomalous signal preserved. On by default with decay/absorption (`--no-scaling-corrections` disables all). The correction-surface cross-validation now scores a sigma-independent R-meas-like agreement instead of a studentized chi^2, so a surface can no longer pass by reshaping the sigmas (which also hardens the absorption surface against a mis-indexed-data regression). ### 1.0.0-rc.159 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. diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index cd2e0ac7..11503447 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -89,8 +89,8 @@ resolution) sharpens both the space-group search and the error model. To tune th `--two-pass-rotation=100` (or `-R100` — the first-pass image count); to force the sweep to be treated as independent stills use `--force-still`. -After the per-frame scale-fulls step, rotation scaling applies two **correction surfaces**, **on by -default** (`--no-scaling-corrections` disables both): +After the per-frame scale-fulls step, rotation scaling applies three **correction surfaces**, **on by +default** (`--no-scaling-corrections` disables all): - **Decay** — a global Debye–Waller relative-*B* over the run, for the radiation damage that weakens later frames more at high resolution (a resolution×time systematic the resolution-flat per-frame @@ -100,11 +100,19 @@ default** (`--no-scaling-corrections` disables both): low photon energy. Its benefit shows up most on model-based metrics: a smooth absorption error largely cancels among symmetry mates (little effect on the error model / ISa) but still biases the intensities, so it measurably lowers *R*free. +- **Modulation** — a smooth multiplicative factor over the position where a reflection lands on the + detector (a flat-field: detector-response and geometric systematics that vary across the detector + plane). Symmetry-equivalents of one reflection land at different detector positions as the crystal + rotates, which over-determines the surface. Because it lives in the detector frame (not the + rotation) the same correction concept applies to stills. This is the largest of the three on + JUNGFRAU data — it lowers *R*meas by several to tens of percent on datasets that carry a + detector systematic, while holding or improving CC1/2 and the anomalous signal. -Both are **cross-validated** — fitted on even-numbered frames and kept only if they improve the -held-out odd-frame symmetry-equivalent agreement by a clear margin (and vice versa) — so where the -systematic is absent they are a no-op rather than a source of added noise; that is why they are safe -to leave on. They run on the GPU when one is present, at negligible cost. +All three are **cross-validated** — fitted on even-numbered frames and kept only if they improve the +held-out odd-frame symmetry-equivalent agreement by a clear margin (and vice versa). The agreement is +scored as a σ-independent, *R*meas-like fractional deviation, so a surface can never pass +cross-validation by merely reshaping the sigmas; where the systematic is absent the surface is a no-op +rather than a source of added noise, which is why they are safe to leave on. ### Still / serial data diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index af6db35d..813782e4 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -171,6 +171,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, // so a no-op when their systematic is absent). Decoupled from the stills-only -B / RefineB flag. refine_decay_b = s.GetCorrectionSurfaces(); absorption_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; + modulation_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; scaling_iter = std::max(1, scaling_iterations); resolution_cutoff_method = s.GetResolutionCutoff(); resolution_cc_target = s.GetResolutionCCTarget(); @@ -661,7 +662,6 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) { return; constexpr int NB = 8; - const int ncell = NB * NB; std::vector cell(fulls.size(), -1); for (size_t i = 0; i < fulls.size(); ++i) { const Obs &o = fulls[i]; @@ -673,14 +673,46 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) { const int iy = std::clamp(static_cast((u.y + 1.0f) * 0.5f * NB), 0, NB - 1); cell[i] = ix * NB + iy; } + ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Absorption (goniometer-frame 8x8)"); +} +void RotationScaleMerge::RefineModulation(int n_iter, int n_groups) { + // Detector-plane modulation (flat-field): a smooth multiplicative factor over the position where a + // reflection lands on the detector (predicted px, py). It absorbs detector-response and geometric + // systematics that vary across the detector plane; symmetry-equivalents of one reflection land at + // different positions as the crystal rotates, over-determining the surface. Because it is a property of + // the detector + beam (not the rotation), the same correction concept transfers to stills. + if (fulls.empty()) + return; + float pxmin = std::numeric_limits::infinity(), pxmax = -pxmin, pymin = pxmin, pymax = -pxmin; + for (const Obs &o : fulls) { + if (!std::isfinite(o.px) || !std::isfinite(o.py)) continue; + pxmin = std::min(pxmin, o.px); pxmax = std::max(pxmax, o.px); + pymin = std::min(pymin, o.py); pymax = std::max(pymax, o.py); + } + if (!(pxmax > pxmin) || !(pymax > pymin)) + return; + constexpr int NB = 16; + const float sx = NB / (pxmax - pxmin), sy = NB / (pymax - pymin); + std::vector cell(fulls.size(), -1); + for (size_t i = 0; i < fulls.size(); ++i) { + const Obs &o = fulls[i]; + if (!std::isfinite(o.px) || !std::isfinite(o.py)) continue; + const int ix = std::clamp(static_cast((o.px - pxmin) * sx), 0, NB - 1); + const int iy = std::clamp(static_cast((o.py - pymin) * sy), 0, NB - 1); + cell[i] = ix * NB + iy; + } + ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Modulation (detector-frame 16x16)"); +} + +void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int ncell, int n_iter, + int n_groups, const char *name) { auto usable = [&](const Obs &o) { return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && o.partiality >= min_partiality; }; - // Fit the per-cell absorption 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 accumulated - // per-cell factor A[cell]. + // 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. auto fit_surface = [&](int parity) -> std::vector { std::vector A(ncell, 1.0); for (int it = 0; it < n_iter; ++it) { @@ -725,39 +757,40 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) { const double sc = static_cast(o.sigma) * o.corr * a, w = 1.0 / (sc * sc); sw[o.group] += w; swI[o.group] += w * Is; } - double sum = 0.0; size_t n = 0; + // 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. + double num = 0.0, den = 0.0; for (size_t i = 0; i < fulls.size(); ++i) { const Obs &o = fulls[i]; if (!usable(o) || cell[i] < 0 || (o.frame & 1) != parity || sw[o.group] <= 0.0) continue; const double a = A[cell[i]], Is = static_cast(o.I) * o.corr * a; - const double sc = static_cast(o.sigma) * o.corr * a; const double Iref = swI[o.group] / sw[o.group]; - if (!std::isfinite(Iref) || !(sc > 0.0)) continue; - const double dd = (Is - Iref) / sc; - sum += dd * dd; ++n; + if (!std::isfinite(Iref) || Iref <= 0.0) continue; + num += std::abs(Is - Iref); den += Iref; } - return n > 0 ? sum / static_cast(n) : 0.0; + return den > 0.0 ? num / den : 0.0; }; // Cross-validate: fit the surface on even frames and score the held-out odd equivalents (and vice - // versa). A real absorption 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 or act on marginal noise. + // 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("Absorption correction: not cross-validated (held-out gain {:.1f}%, skipped)", - 100.0 * gain / std::max(base, 1e-30)); + logger.Info("{} correction: not cross-validated (held-out gain {:.1f}%, skipped)", + name, 100.0 * gain / std::max(base, 1e-30)); return; } const std::vector A = fit_surface(-1); for (size_t i = 0; i < fulls.size(); ++i) if (cell[i] >= 0) fulls[i].corr = static_cast(fulls[i].corr * A[cell[i]]); - logger.Info("Absorption correction: goniometer-frame {}x{} surface (cross-validated, held-out gain {:.1f}%)", - NB, NB, 100.0 * gain / std::max(base, 1e-30)); + logger.Info("{} correction: cross-validated, held-out gain {:.1f}%", + name, 100.0 * gain / std::max(base, 1e-30)); } void RotationScaleMerge::ComputeSmoothGWindow(const std::vector &g, int window, @@ -909,7 +942,7 @@ void RotationScaleMerge::Combine() { if (dump != nullptr) *dump << full.h << ' ' << full.k << ' ' << full.l << ' ' << full.I << ' ' << full.sigma << ' ' << d << ' ' << n_frames_event << ' ' << sum_partiality << ' ' - << static_cast(peak_frame) << '\n'; + << static_cast(peak_frame) << ' ' << full.px << ' ' << full.py << '\n'; } return ev.size(); }; @@ -922,7 +955,9 @@ void RotationScaleMerge::Combine() { std::ofstream dump; if (!observation_dump_path.empty()) { dump.open(observation_dump_path); - dump << "# h k l I sigma d n_frames captured_fraction\n"; + // Columns are enough to rescale/merge the combined fulls in an external program: peak_frame + // is a per-full batch (ZD) and px/py the predicted detector position (XD/YD). + dump << "# h k l I sigma d n_frames captured_fraction peak_frame px py\n"; } fulls.reserve(n_run); for (int r = 0; r < n_run; ++r) @@ -1600,11 +1635,14 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, // the merged reference. Cheap host loops over the downloaded fulls; applied only on the final // in-symmetry merge, never the P1 space-group search pass (corrections there add risk and can // perturb the symmetry determination). --- - [[maybe_unused]] const bool corrections = !for_search && (refine_decay_b || absorption_iter > 0); + [[maybe_unused]] const bool corrections = + !for_search && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0); if (!for_search && refine_decay_b) RefineDecay(n_groups); if (!for_search && absorption_iter > 0) RefineAbsorption(absorption_iter, n_groups); + if (!for_search && modulation_iter > 0) + RefineModulation(modulation_iter, n_groups); #ifdef JFJOCH_USE_CUDA // The corrections mutate the host fulls' corr; when the merge runs on the resident (GPU) fulls, push // the corrected corr back to the device so the merge reads it. diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 8bac7791..f92da2bf 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -106,6 +106,7 @@ private: bool scale_fulls = true; bool refine_decay_b = false; // per-time-block Debye-Waller decay correction (radiation damage) int absorption_iter = 0; // >0: fit a goniometer-frame absorption surface over this many iterations + int modulation_iter = 0; // >0: fit a detector-plane modulation (flat-field) surface, this many iterations double mosaicity_deg = 0.1; float ice_half_width_q = 0.0f; // Automatic high-resolution cutoff for the written reflections + reported shells (post-merge; the @@ -198,6 +199,16 @@ private: // the goniometer frame (path-length / absorption; negligible at hard X-rays, matters at low energy). void RefineDecay(int n_groups); void RefineAbsorption(int n_iter, int n_groups); + // Detector-plane modulation (flat-field): the same cross-validated surface fit as absorption, but the + // cell is the predicted detector position (px, py) instead of the goniometer-frame direction. Corrects + // detector-response / geometric systematics that vary with where a reflection lands; because it lives + // 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. + void ApplyCellSurface(const std::vector &cell, int ncell, int n_iter, int n_groups, + const char *name); // 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.