From c5cfdcef7ea537a988d90d8b3921ca03d8f0367e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 18 Jul 2026 22:08:09 +0200 Subject: [PATCH] rugnux: per-batch relative-B correction + radiation-damage B monitor Rotation scaling modelled dose-dependent scattering power with a single global relative-B decay slope only (RefineDecay); per-frame refinement was scale-only. Add a per-batch relative Debye-Waller B, both as an opt-in correction and as an always-on measurement. - FitRelativeBCurve / SolveCurvatureSmoothedB: per rotation-range batch, fit the s^2 slope of ln(Iref/Iobs) that the resolution-flat per-frame G leaves, regularised by a second-difference (curvature) penalty so a genuine smooth relative-B is kept while resolution-correlated noise (which an unconstrained per-batch fit would chase into an oscillating curve) is suppressed. - Correction (--relative-b[=deg], default off): applied after RefineDecay, gated by a physical peak-to-peak floor and cross-validated by ASU-group parity (a per-batch parameter owns whole frames, so it is scored on held-out equivalents, not held-out frames). On a battery re-scaling identical stored reflections it is a no-op on well-behaved data and a small gain on radiation-damaged rotation data; never a regression. - Radiation-damage monitor (always on for rotation, report-only): measure the relative-B change from the pristine low-dose start of the run to its end (a whole-run-mean reference inverts on damaged data) plus the per-batch curve, and report them on the command line (in the radiation-damage report and the merge statistics) and in the mmCIF (_reflns.jfjoch_radiation_damage_ relative_B and a per-batch loop). A large magnitude flags radiation damage, complementing the existing per-image CC / mosaicity read. - --scale now reproduces the de-novo rotation merge (it previously omitted scale-fulls and smooth-G, silently using a weaker model), and honours --relative-b, so an offline re-scale matches the full pipeline. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/ScalingSettings.cpp | 11 + common/ScalingSettings.h | 9 + image_analysis/WriteReflections.cpp | 20 ++ image_analysis/scale_merge/Merge.cpp | 4 + image_analysis/scale_merge/Merge.h | 9 + .../scale_merge/RotationScaleMerge.cpp | 245 +++++++++++++++++- .../scale_merge/RotationScaleMerge.h | 29 +++ rugnux/Rugnux.cpp | 13 + rugnux/rugnux_cli.cpp | 22 ++ 9 files changed, 361 insertions(+), 1 deletion(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index a86492ae..01477cb9 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -171,6 +171,17 @@ double ScalingSettings::GetSmoothGDegrees() const { return smooth_g_deg; } +ScalingSettings &ScalingSettings::RelativeBDegrees(double input) { + if (input < 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Relative-B batch width must be non-negative"); + relative_b_deg = input; + return *this; +} + +double ScalingSettings::GetRelativeBDegrees() const { + return relative_b_deg; +} + double ScalingSettings::GetMinPartiality() const { return min_partiality; } diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index d4c9ec6e..74105de7 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -57,6 +57,13 @@ class ScalingSettings { // the smoothing physical (independent of frame slicing). 0 = off. A no-op without rot3d. double smooth_g_deg = 0.0; + // Per-batch relative-B on the rot3d fulls (beyond the single global decay slope): bin frames into + // rotation-range batches of this width in degrees and refine one relative Debye-Waller B per batch, so + // NON-monotonic changes in scattering power across the run (absorption path, crystal slippage, dose + // bursts) are corrected the resolution-flat per-frame G and the single decay slope both miss. Cross- + // validated + zero-mean-anchored, so a no-op when absent. 0 = off. A no-op without rot3d. + double relative_b_deg = 0.0; + double rfree_fraction = 0.05; // Automatic high-resolution cutoff for the written reflections + reported shells (not the scaling @@ -86,6 +93,7 @@ public: ScalingSettings& CorrectionSurfaces(bool input); ScalingSettings& StillsModulation(bool input); ScalingSettings& SmoothGDegrees(double input); + ScalingSettings& RelativeBDegrees(double input); ScalingSettings& RfreeFraction(double input); ScalingSettings& ScalingRegularize(bool input); @@ -124,6 +132,7 @@ public: [[nodiscard]] bool GetCorrectionSurfaces() const; [[nodiscard]] bool GetStillsModulation() const; [[nodiscard]] double GetSmoothGDegrees() const; + [[nodiscard]] double GetRelativeBDegrees() const; [[nodiscard]] double GetRfreeFraction() const; [[nodiscard]] bool GetScalingRegularize() const; diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 80017e97..fffb4ab1 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -170,8 +170,28 @@ void WriteMmcifReflections(const std::vector &reflections, if (twinning.moment_reflections > 0) out << "_reflns.jfjoch_second_moment_I " << Fmt(twinning.second_moment, 3) << " # /^2 (untwinned 2.00, perfect twin 1.50)\n"; + // Radiation-damage monitor (rotation): the relative Debye-Waller B change from the first to the last + // frame (A^2). A large magnitude flags a dose-dependent resolution-scale change = radiation damage; + // positive is the typical direction (high-resolution intensity fades with dose). No standard mmCIF item. + if (std::isfinite(statistics.radiation_damage_delta_b)) + out << "_reflns.jfjoch_radiation_damage_relative_B " << Fmt(statistics.radiation_damage_delta_b, 2) + << " # relative-B first->last over the run (A^2); + = high-res fades with dose\n"; out << "#\n"; + // Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B + // per rotation-range batch, measured before any correction. rotation_start_deg = id * batch_deg. + if (!statistics.radiation_damage_b_batch.empty()) { + out << "loop_\n"; + out << "_jfjoch_radiation_damage_batch.id\n"; + out << "_jfjoch_radiation_damage_batch.rotation_start_deg\n"; + out << "_jfjoch_radiation_damage_batch.relative_B\n"; + for (size_t i = 0; i < statistics.radiation_damage_b_batch.size(); ++i) + out << " " << (i + 1) << " " + << Fmt(static_cast(i) * statistics.radiation_damage_batch_deg, 1) << " " + << Fmt(statistics.radiation_damage_b_batch[i], 2) << "\n"; + out << "#\n"; + } + out << "loop_\n"; out << "_reflns_shell.d_res_high\n"; out << "_reflns_shell.d_res_low\n"; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 518fee26..b67484ef 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -923,6 +923,10 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { if (std::isfinite(in.wilson_b) && in.wilson_b > 0.0) output << fmt::format(" Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f})", in.wilson_b, in.wilson_b_correlation) << std::endl; + if (std::isfinite(in.radiation_damage_delta_b)) + output << fmt::format(" Radiation damage: relative B-factor change over run = {:+.2f} A^2 " + "(first->last, {} batches)", + in.radiation_damage_delta_b, in.radiation_damage_b_batch.size()) << std::endl; output << std::endl; return output; } diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 2636a1f7..441252c9 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -62,6 +62,15 @@ struct MergeStatistics { // Diagnostic only; not used in scaling. NaN when not determined. double wilson_b = NAN; double wilson_b_correlation = NAN; + + // Radiation-damage monitor (rotation only): the relative Debye-Waller B change measured from the first + // to the last frame of the run (A^2; positive = high-resolution intensity fades with dose = damage) and + // the per-batch relative-B curve it was derived from. Measured before any decay/relative-B correction is + // applied, so it reports how much radiation damage was present. Diagnostic; NaN / empty for stills or + // when not determined. batch_deg is the rotation width per batch of the curve. + double radiation_damage_delta_b = NAN; + std::vector radiation_damage_b_batch; + double radiation_damage_batch_deg = 0.0; }; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 80d1b6d0..27012c7a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -173,6 +173,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, refine_decay_b = s.GetCorrectionSurfaces(); absorption_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; modulation_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; + relative_b_deg = s.GetRelativeBDegrees(); // opt-in per-batch relative-B (0 = off) scaling_iter = std::max(1, scaling_iterations); resolution_cutoff_method = s.GetResolutionCutoff(); resolution_cc_target = s.GetResolutionCCTarget(); @@ -650,6 +651,236 @@ void RotationScaleMerge::RefineDecay(int n_groups) { 0.5 * slope * n_frames, 0.5 * slope); } +std::vector RotationScaleMerge::SolveCurvatureSmoothedB(const std::vector &num, + const std::vector &den) const { + // Solve for a SMOOTH per-batch relative-B from the per-batch normal-equation data (num_c = sum w s^2 y, + // den_c = sum w s^4). A per-batch regression alone (b_c = num_c/den_c) over-fits: a genuine relative-B is + // a low-curvature function of dose / orientation, but an unconstrained fit chases resolution-correlated + // noise into a rough, oscillating curve. So minimise data-fidelity PLUS a second-difference (curvature) + // penalty sum_c den_c (b_c - num_c/den_c)^2 + mu * sum_c (b_{c-1} - 2 b_c + b_{c+1})^2, which leaves a + // straight RAMP untouched but suppresses oscillation, solving the small SPD normal system A b = num, + // A = diag(den + lambda) + mu * D^T D, by Gauss-Seidel (n_batch <= a few tens). Not anchored - callers set + // the gauge (zero-mean for the correction, low-dose reference for the monitor). + constexpr double RELATIVE_B_MAX = 50.0; // A^2, clamp a single batch's relative-B (generous guard) + constexpr double RELATIVE_B_CURVATURE_K = 8.0; // curvature penalty weight (x the median batch leverage) + const int n_batch = static_cast(num.size()); + std::vector dsorted(den); + std::sort(dsorted.begin(), dsorted.end()); + const double med_den = dsorted.empty() ? 0.0 : dsorted[dsorted.size() / 2]; + const double lambda = 0.1 * med_den; // ridge ~ 10% of the median batch leverage + const double mu = RELATIVE_B_CURVATURE_K * med_den; // curvature penalty + std::vector A(static_cast(n_batch) * n_batch, 0.0); + auto Aat = [&](int i, int j) -> double & { return A[static_cast(i) * n_batch + j]; }; + for (int c = 0; c < n_batch; ++c) Aat(c, c) = den[c] + lambda; + constexpr int st[3] = {-1, 0, 1}; + constexpr double sv[3] = {1.0, -2.0, 1.0}; // second-difference stencil + for (int c = 1; c < n_batch - 1; ++c) + for (int a = 0; a < 3; ++a) + for (int b2 = 0; b2 < 3; ++b2) + Aat(c + st[a], c + st[b2]) += mu * sv[a] * sv[b2]; + std::vector b(n_batch, 0.0); + for (int it = 0; it < 100; ++it) + for (int c = 0; c < n_batch; ++c) { + if (!(Aat(c, c) > 0.0)) continue; + double s = num[c]; + for (int j = 0; j < n_batch; ++j) + if (j != c) s -= Aat(c, j) * b[j]; + b[c] = std::clamp(s / Aat(c, c), -RELATIVE_B_MAX, RELATIVE_B_MAX); + } + return b; +} + +std::vector RotationScaleMerge::FitRelativeBCurve(int n_groups, int n_batch, int frames_per_batch, + int gparity) const { + // Fit one relative-B per batch over the ASU-group subset {group&1 == gparity} (gparity < 0 = all fulls), + // against an inverse-variance reference built from that same subset (no leakage). Per batch the model is + // y = ln(I_ref / I_obs) = b * s^2 through the origin (the resolution-flat intercept is the per-frame G, + // already applied), so the data term is b_c = sum(w s^2 y) / sum(w s^2 s^2); SolveCurvatureSmoothedB then + // smooths it. Finally subtract the leverage-weighted mean: a constant B is a global Wilson-B degenerate + // with overall scale, and cannot change equivalent agreement (within one ASU group every equivalent sits + // at the same s^2, so only b VARYING between batches moves the merge). Drives the per-batch correction. + auto s2_of = [](float d) { return d > 0.0f ? 1.0 / (4.0 * static_cast(d) * d) : 0.0; }; + auto usable = [&](const Obs &o) { + return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && (o.d > 0.0f) + && o.partiality >= min_partiality; + }; + auto batch_of = [&](const Obs &o) { return std::min(n_batch - 1, o.frame / frames_per_batch); }; + + std::vector sw(n_groups, 0.0), swI(n_groups, 0.0); + for (const auto &o : fulls) { + if (!usable(o) || (gparity >= 0 && (o.group & 1) != gparity)) continue; + const double sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); + sw[o.group] += w; swI[o.group] += w * static_cast(o.I) * o.corr; + } + std::vector num(n_batch, 0.0), den(n_batch, 0.0); + for (const auto &o : fulls) { + if (!usable(o) || (gparity >= 0 && (o.group & 1) != gparity) || sw[o.group] <= 0.0) continue; + const double Iref = swI[o.group] / sw[o.group]; + const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; + if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue; + const double w = (Is / sc) * (Is / sc); + const double s2 = s2_of(o.d), y = std::log(Iref / Is); + num[batch_of(o)] += w * s2 * y; den[batch_of(o)] += w * s2 * s2; + } + std::vector b = SolveCurvatureSmoothedB(num, den); + double dw = 0.0, dbw = 0.0; + for (int c = 0; c < n_batch; ++c) { dw += den[c]; dbw += den[c] * b[c]; } + const double bmean = dw > 0.0 ? dbw / dw : 0.0; + for (auto &v : b) v -= bmean; + return b; +} + +void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { + // Radiation-damage MONITOR (report-only): measure the per-batch relative Debye-Waller B on the scaled + // fulls BEFORE any correction, so it captures how the high-resolution scattering power fades with + // accumulated dose. Crucially the reference is the PRISTINE LOW-DOSE crystal (first ~10% of frames), not + // the whole-run mean: a group measured both early and late then contributes ln(I_early / I_late) vs s^2, + // whose slope is the ADDED B - an absolute, correctly-signed damage read (positive = high-res fades). + // A whole-run-mean reference inverts on heavily damaged data (the surviving late high-res reflections + // look strong against a mean dominated by the frames where they had already faded). Stores the first->last + // relative-B change (the headline "B-factor difference") and the per-batch curve; never touches corr. + rad_damage_delta_b = std::numeric_limits::quiet_NaN(); + rad_damage_b_batch.clear(); + rad_damage_batch_deg = 0.0; + if (fulls.empty() || n_frames <= 0) + return; + const auto gon = x.GetGoniometer(); + const double osc_deg = gon ? std::fabs(gon->GetIncrement_deg()) : 0.0; + if (!(osc_deg > 1e-6)) + return; + constexpr double MONITOR_BATCH_DEG = 10.0; // rotation width per batch for the damage curve + const int frames_per_batch = std::max(1, static_cast(std::lround(MONITOR_BATCH_DEG / osc_deg))); + const int n_batch = std::max(1, (n_frames + frames_per_batch - 1) / frames_per_batch); + if (n_batch < 2) // a single batch is a global Wilson-B (degenerate with scale) - nothing to monitor + return; + auto s2_of = [](float d) { return d > 0.0f ? 1.0 / (4.0 * static_cast(d) * d) : 0.0; }; + auto usable = [&](const Obs &o) { + return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && (o.d > 0.0f) + && o.partiality >= min_partiality; + }; + auto batch_of = [&](const Obs &o) { return std::min(n_batch - 1, o.frame / frames_per_batch); }; + + // Low-dose reference: inverse-variance group means over the first ~10% of frames (>= 1 batch). + const int ref_frames = std::max(frames_per_batch, n_frames / 10); + std::vector sw0(n_groups, 0.0), swI0(n_groups, 0.0); + for (const auto &o : fulls) { + if (!usable(o) || o.frame >= ref_frames) continue; + const double sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); + sw0[o.group] += w; swI0[o.group] += w * static_cast(o.I) * o.corr; + } + // Per batch: slope of ln(I_ref0 / I_obs) vs s^2, over groups also seen in the low-dose reference. + std::vector num(n_batch, 0.0), den(n_batch, 0.0); + for (const auto &o : fulls) { + if (!usable(o) || sw0[o.group] <= 0.0) continue; + const double Iref = swI0[o.group] / sw0[o.group]; + const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; + if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue; + const double w = (Is / sc) * (Is / sc); + const double s2 = s2_of(o.d), y = std::log(Iref / Is); + num[batch_of(o)] += w * s2 * y; den[batch_of(o)] += w * s2 * s2; + } + std::vector b = SolveCurvatureSmoothedB(num, den); + // Anchor so the reference batches read ~0 (b = added B relative to the low-dose start). + const int n_ref_batch = std::max(1, ref_frames / frames_per_batch); + double ref_mean = 0.0; + for (int c = 0; c < n_ref_batch && c < n_batch; ++c) ref_mean += b[c]; + ref_mean /= std::min(n_ref_batch, n_batch); + for (auto &v : b) v -= ref_mean; + // Headline: the first->last relative-B change from a least-squares linear trend of the (smoothed) curve - + // robust to a single noisy end batch, unlike raw endpoints. Positive = high-res fades with dose = damage. + double sx = 0, sy = 0, sxx = 0, sxy = 0; + for (int c = 0; c < n_batch; ++c) { sx += c; sy += b[c]; sxx += double(c) * c; sxy += double(c) * b[c]; } + const double det = n_batch * sxx - sx * sx; + const double slope = det > 0.0 ? (n_batch * sxy - sx * sy) / det : 0.0; + rad_damage_delta_b = slope * (n_batch - 1); + rad_damage_batch_deg = MONITOR_BATCH_DEG; + rad_damage_b_batch.resize(n_batch); + for (int c = 0; c < n_batch; ++c) rad_damage_b_batch[c] = static_cast(b[c]); +} + +void RotationScaleMerge::RefineRelativeB(int n_groups) { + // RefineDecay removes the AVERAGE radiation-damage falloff as a single global relative-B slope, but the + // relative scattering power drifts NON-monotonically across a run (absorption path as the crystal + // rotates, slippage, dose bursts, a partly-shadowed wedge). XDS/AIMLESS model this with a per-BATCH + // relative-B on top of the per-batch scale. Here, on the fulls after RefineDecay: fit the smoothed + // per-batch relative-B (FitRelativeBCurve) and apply exp(b_batch * s^2). Same anti-overfit posture as + // RefineDecay - a physical spread floor and a cross-validated held-out gain - but the CV split is by ASU + // GROUP parity, not frame parity: a per-batch parameter owns whole frames, so it cannot be scored on a + // held-out FRAME; splitting the equivalents instead tests whether a batch's B generalises to reflections + // it was not fit on, which is exactly the over-fit risk of a many-parameter model. + if (fulls.empty() || n_frames <= 0 || relative_b_deg <= 0.0) + return; + const auto gon = x.GetGoniometer(); + const double osc_deg = gon ? std::fabs(gon->GetIncrement_deg()) : 0.0; + if (!(osc_deg > 1e-6)) + return; + const int frames_per_batch = std::max(1, static_cast(std::lround(relative_b_deg / osc_deg))); + const int n_batch = std::max(1, (n_frames + frames_per_batch - 1) / frames_per_batch); + if (n_batch < 2) // a single batch is a global Wilson-B (degenerate with overall scale) - nothing relative + return; + auto s2_of = [](float d) { return d > 0.0f ? 1.0 / (4.0 * static_cast(d) * d) : 0.0; }; + auto usable = [&](const Obs &o) { + return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && (o.d > 0.0f) + && o.partiality >= min_partiality; + }; + auto batch_of = [&](const Obs &o) { return std::min(n_batch - 1, o.frame / frames_per_batch); }; + auto b_factor = [&](const Obs &o, const std::vector &b) { + return std::exp(b[batch_of(o)] * s2_of(o.d)); + }; + // Held-out disagreement over the group subset {group&1 == gparity} with `b` applied: the same + // sigma-INDEPENDENT Rmeas-like metric as RefineDecay (sum|Is-Iref|/sum|Iref|), so a relative-B cannot + // "win" by reshaping sigma. Reference built from the scored subset (no leakage). + auto subset_disagreement = [&](int gparity, const std::vector &b) -> double { + std::vector sw(n_groups, 0.0), swI(n_groups, 0.0); + for (const auto &o : fulls) { + if (!usable(o) || (o.group & 1) != gparity) continue; + const double f = b_factor(o, b); + const double Is = static_cast(o.I) * o.corr * f, sc = static_cast(o.sigma) * o.corr * f; + const double w = 1.0 / (sc * sc); + sw[o.group] += w; swI[o.group] += w * Is; + } + double num = 0.0, den = 0.0; + for (const auto &o : fulls) { + if (!usable(o) || (o.group & 1) != gparity || sw[o.group] <= 0.0) continue; + const double Iref = swI[o.group] / sw[o.group]; + const double Is = static_cast(o.I) * o.corr * b_factor(o, b); + if (!std::isfinite(Iref) || Iref <= 0.0) continue; + num += std::abs(Is - Iref); den += Iref; + } + return den > 0.0 ? num / den : 0.0; + }; + + // A relative-B whose peak-to-peak is below a couple A^2 is negligible; "correcting" it only spreads + // equivalents (same s^2, different batch) into extra scatter - net harmful, exactly like RefineDecay's + // sub-floor slope. (RELATIVE_B_MIN_SPREAD ~ DECAY_MIN_DELTA_B.) + constexpr double RELATIVE_B_MIN_SPREAD = 2.0; // A^2, minimum peak-to-peak relative-B to engage + const std::vector b_all = FitRelativeBCurve(n_groups, n_batch, frames_per_batch, -1); + const double bmin = *std::min_element(b_all.begin(), b_all.end()); + const double bmax = *std::max_element(b_all.begin(), b_all.end()); + if (bmax - bmin < RELATIVE_B_MIN_SPREAD) { + logger.Info("Relative-B: negligible variation ({} batches, peak-to-peak {:.2f} A^2 < {:.1f}, skipped)", + n_batch, bmax - bmin, RELATIVE_B_MIN_SPREAD); + return; + } + // Cross-validate: fit on even ASU groups, score the held-out odd equivalents (and vice-versa). Real + // per-batch structure generalises across the equivalent split; over-fit per-batch noise does not. + const std::vector zero(n_batch, 0.0); + const double base = subset_disagreement(0, zero) + subset_disagreement(1, zero); + const std::vector b_even = FitRelativeBCurve(n_groups, n_batch, frames_per_batch, 0); + const std::vector b_odd = FitRelativeBCurve(n_groups, n_batch, frames_per_batch, 1); + const double gain = base - (subset_disagreement(1, b_even) + subset_disagreement(0, b_odd)); + if (!(gain > CV_MIN_RELATIVE_GAIN * base)) { + logger.Info("Relative-B: not cross-validated ({} batches, peak-to-peak {:.2f} A^2, held-out gain {:.1f}%, skipped)", + n_batch, bmax - bmin, 100.0 * gain / std::max(base, 1e-30)); + return; + } + for (auto &o : fulls) + if (usable(o)) + o.corr = static_cast(o.corr * b_factor(o, b_all)); + logger.Info("Relative-B: {} batches of ~{:.1f} deg, peak-to-peak {:.2f} A^2, held-out gain {:.1f}% (cross-validated)", + n_batch, relative_b_deg, bmax - bmin, 100.0 * gain / std::max(base, 1e-30)); +} + void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) { // Absorption / path-length: a smooth multiplicative factor over the diffracted-beam direction in the // goniometer (crystal) frame. Each full's lab diffracted direction (from its predicted detector @@ -1567,6 +1798,12 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool logger.Info("Anomalous signal SigAno = {:.2f} ({} acentric pairs)", overall.abs_diff_over_sigma_anomalous, sig_n); + // Radiation-damage monitor (measured before any correction by MeasureRadiationDamageB): carry the + // first->last relative-B change and the per-batch curve into the reported statistics / mmCIF. + out.radiation_damage_delta_b = rad_damage_delta_b; + out.radiation_damage_b_batch = rad_damage_b_batch; + out.radiation_damage_batch_deg = rad_damage_batch_deg; + // Attach the per-reflection anomalous split so the writer can emit I(+)/I(-) by default (each merged // reflection maps to its Friedel-ASU key; in an anomalous merge both mates map to the same key). if (!anom_export.empty()) { @@ -1757,9 +1994,15 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, // 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 || modulation_iter > 0); + !for_search && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0); + // Radiation-damage monitor: measure the per-batch relative-B on the scaled fulls BEFORE any correction + // (report-only; captures the full damage signature, not a residual). Skipped on the P1 search pass. + if (!for_search) + MeasureRadiationDamageB(n_groups); if (!for_search && refine_decay_b) RefineDecay(n_groups); + if (!for_search && relative_b_deg > 0.0) + RefineRelativeB(n_groups); // per-batch relative-B on top of the single decay slope if (!for_search && absorption_iter > 0) RefineAbsorption(absorption_iter, n_groups); if (!for_search && modulation_iter > 0) diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index f92da2bf..1ca5c5a5 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -4,6 +4,7 @@ #pragma once #include +#include #include #include @@ -107,6 +108,7 @@ private: 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 relative_b_deg = 0.0; // >0: fit a per-batch relative-B (batch width in deg); 0 = off 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 @@ -142,6 +144,12 @@ private: // written back for the per-image scaling table. Empty if there is no per-frame mosaicity. std::vector mos_smooth; + // Radiation-damage monitor (measured by MeasureRadiationDamageB on the scaled fulls before any decay + // correction; report-only, copied into the result statistics by MergeAndStats). NaN / empty until set. + double rad_damage_delta_b = std::numeric_limits::quiet_NaN(); // relative-B first->last (A^2) + std::vector rad_damage_b_batch; // per-batch relative-B curve (A^2) + double rad_damage_batch_deg = 0.0; // rotation width per batch (deg) + // Working per-group arrays (sized to the current group count; reused). std::vector group_h, group_k, group_l; @@ -198,6 +206,27 @@ private: // physical total-dB floor). RefineAbsorption fits a smooth factor over the diffracted-beam direction in // the goniometer frame (path-length / absorption; negligible at hard X-rays, matters at low energy). void RefineDecay(int n_groups); + // Solve a smooth per-batch relative-B from the per-batch normal equations (num_c=sum w s^2 y, + // den_c=sum w s^4): data-fidelity + a second-difference (curvature) penalty, by Gauss-Seidel. Returns the + // un-anchored curve; the caller sets the gauge. Shared by the correction and the radiation-damage monitor. + std::vector SolveCurvatureSmoothedB(const std::vector &num, + const std::vector &den) const; + // Fit a smoothed per-batch relative-B curve (A^2 per batch) on the fulls over the ASU-group subset + // {group&1==gparity} (gparity<0 = all): the weighted s^2 slope of ln(Iref/Iobs) per batch against a + // subset-global reference, smoothed and zero-mean-anchored. Drives the per-batch correction. + std::vector FitRelativeBCurve(int n_groups, int n_batch, int frames_per_batch, int gparity) const; + // Radiation-damage MONITOR (report-only): measure the per-batch relative-B on the scaled fulls before + // any decay correction and store the first->last relative-B change + the per-batch curve on this object + // (copied into the result statistics by MergeAndStats, then printed / logged / written to the mmCIF). + void MeasureRadiationDamageB(int n_groups); + // Per-batch relative-B, applied after RefineDecay: the single decay slope removes the average + // radiation-damage falloff, but the relative scattering power drifts NON-monotonically across a run + // (absorption path, crystal slippage, dose bursts). Refine one relative Debye-Waller B per batch + // (FitRelativeBCurve), anchored to zero mean (the constant part is a global Wilson-B, degenerate with + // overall scale). Guarded by a physical peak-to-peak floor and cross-validated by ASU-GROUP parity (a + // per-batch parameter cannot be scored on a held-out FRAME the batch owns; splitting the equivalents + // tests whether a batch's B generalises to reflections it was not fit on). Opt-in (--relative-b). + void RefineRelativeB(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 diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 862898ba..826ba950 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1389,6 +1389,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b cc0, cc9, cc9 < cc0 - 0.05 ? " (falling: radiation damage / crystal decay)" : ""); } + // Relative B-factor across the sweep (measured before any decay correction, against the + // low-dose start): the resolution-dependent complement to the CC/mosaicity read above. A large + // magnitude flags radiation damage; positive is the typical direction (high-res fades with dose). + if (std::isfinite(sm.statistics.radiation_damage_delta_b)) { + os << fmt::format("\n => relative B-factor change over run = {:+.2f} A^2 (first->last){}", + sm.statistics.radiation_damage_delta_b, + std::fabs(sm.statistics.radiation_damage_delta_b) > 5.0 + ? " (significant dose-dependent scaling: radiation damage)" : ""); + os << fmt::format("\n per-batch relative-B (A^2, {:.0f} deg/batch):", + sm.statistics.radiation_damage_batch_deg); + for (float bb : sm.statistics.radiation_damage_b_batch) + os << fmt::format(" {:.1f}", bb); + } logger.Info("{}", os.str()); } } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9ae4e4c2..d8bebcab 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -37,6 +37,9 @@ // Default rot3d per-frame scale-G smoothing range (XDS DELPHI-like), in degrees of rotation. constexpr double SMOOTH_G_DEFAULT_DEG = 5.0; +// Default rot3d per-batch relative-B batch width (bare --relative-b), in degrees of rotation. +constexpr double RELATIVE_B_DEFAULT_DEG = 10.0; + // Default rot3d per-observation merge outlier rejection (sigma from the per-reflection median, XDS-like). // Neutral-or-better across the test crystals; pass --reject-outliers 0 to disable. constexpr double REJECT_OUTLIERS_DEFAULT_NSIGMA = 6.0; @@ -88,6 +91,7 @@ void print_usage() { std::cout << " --no-scale-fulls Disable the rot3d scale-fulls refit (it is on by default for rot3d)" << std::endl; std::cout << " --write-process-h5 Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl; std::cout << " --smooth-g[=deg] rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl; + std::cout << " --relative-b[=deg] rot3d: fit a per-batch relative-B (beyond the single decay slope) over deg-degree batches; cross-validated (default: 10 deg when bare; off otherwise)" << std::endl; std::cout << " --no-scaling-corrections rot3d: disable the (default-on) decay + absorption correction surfaces fitted on the fulls after scale-fulls" << std::endl; std::cout << " --stills-modulation stills: fit a detector-plane modulation (flat-field) surface over the merged reflections (cross-validated; experimental, default off)" << std::endl; std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; @@ -169,6 +173,7 @@ enum { OPT_MIN_CAPTURED_FRACTION, OPT_MOSAICITY, OPT_SMOOTH_G, + OPT_RELATIVE_B, OPT_NO_SCALING_CORRECTIONS, OPT_STILLS_MODULATION, OPT_DETECT_ICE_RINGS, @@ -215,6 +220,7 @@ static option long_options[] = { {"no-scale-fulls", no_argument, nullptr, OPT_NO_SCALE_FULLS}, {"write-process-h5", no_argument, nullptr, OPT_WRITE_PROCESS_H5}, {"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G}, + {"relative-b", optional_argument, nullptr, OPT_RELATIVE_B}, {"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS}, {"stills-modulation", no_argument, nullptr, OPT_STILLS_MODULATION}, {"refine", required_argument, nullptr, 'r'}, @@ -497,6 +503,7 @@ int main(int argc, char **argv) { // Geometry overrides (default: keep the value stored in the input file) std::optional beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, polarization_factor; std::optional smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise + std::optional relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption surfaces bool stills_modulation_flag = false; // --stills-modulation: detector-plane flat-field surface for stills bool anomalous_mode = false; @@ -780,6 +787,9 @@ int main(int argc, char **argv) { case OPT_SMOOTH_G: smooth_g_deg_arg = optarg ? parse_double_arg(optarg, "--smooth-g", logger) : SMOOTH_G_DEFAULT_DEG; break; + case OPT_RELATIVE_B: + relative_b_deg_arg = optarg ? parse_double_arg(optarg, "--relative-b", logger) : RELATIVE_B_DEFAULT_DEG; + break; case OPT_STILLS_MODULATION: stills_modulation_flag = true; break; @@ -1068,6 +1078,17 @@ int main(int argc, char **argv) { scaling_settings.OutlierRejectNsigma( outlier_reject_nsigma.value_or( (experiment.GetGoniometer().has_value() && !force_still) ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0)); + // Mirror the de-novo rotation scaling defaults so --scale reproduces the full-pipeline merge on the + // same reflections (scale-fulls + smooth-G, and the opt-in per-batch relative-B). Without these the + // offline re-scale silently used a weaker model than the pipeline that wrote the _process.h5. + { + const bool rot = experiment.GetGoniometer().has_value() && !force_still; + scaling_settings.ScaleFulls(scale_fulls_arg.value_or(rot)); + scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(rot ? SMOOTH_G_DEFAULT_DEG : 0.0)); + scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); + if (no_scaling_corrections) + scaling_settings.CorrectionSurfaces(false); + } experiment.ImportScalingSettings(scaling_settings); if (!experiment.GetUnitCell()) { @@ -1362,6 +1383,7 @@ int main(int argc, char **argv) { ScalingSettings scaling_settings; scaling_settings.ScaleFulls(scale_fulls); scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(rotation_indexing ? SMOOTH_G_DEFAULT_DEG : 0.0)); + scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off if (no_scaling_corrections) scaling_settings.CorrectionSurfaces(false); scaling_settings.StillsModulation(stills_modulation_flag);