From d38100e2152cb0e6307c7a341ef1ddc511b82bb0 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 18:35:47 +0100 Subject: [PATCH] RotationScaleMerge: the decay over sweep segments is a ramp along the dose order or one B per sweep, decided on held-out frames Two hypotheses of how damage accrued over several sweeps: a B rising linearly along the dose order the segments' relative B factors put them in, or one B per sweep - a step at each boundary, flat within a sweep. Both are regressions of ln(I_ref / I_obs) on s^2 within each frame (each frame has its own scale, which took up its part of the B), at the weights the single-sweep slope fit uses. Each is fitted on the even frames and scored by its weighted squared residual on the odd ones, against their own reference, and the other way round; the one that predicts the held-out frames better is applied where it beats no decay by the margin and physical floor a correction has to clear on one sweep. The per-sweep B is anchored at the frame-weighted mean, as the ramp is at the dose centre. The residual is read at weights fixed by the observations, so a correction cannot win by reshaping sigma; the sigma-free disagreement the single-sweep decay is gated on favoured the smaller correction (on three sweeps at B 0, 4 and 8 the right per-sweep B scored worse than a quarter of it), so it is not used for this choice. Result reports the choice (segment_decay_model). Without segments nothing changes. Synthetic test pairs: three sweeps each flat at its own B (8, 0, 4 A^2): held-out residual 33.1% of no decay for the ramp, 22.4% for one B per sweep => per-sweep B. Two sweeps of one steady ramp (6 -> 12 and 0 -> 6, given in that order): ramp 30.7%, per-sweep 47.3% => dose-linear, dB 10.2 A^2. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- .../scale_merge/RotationScaleMerge.cpp | 152 ++++++++++++++++-- .../scale_merge/RotationScaleMerge.h | 7 + tests/RotationScaleMergeSegmentsTest.cpp | 23 ++- 3 files changed, 170 insertions(+), 12 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 694ca83d4..2848ca1c4 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -2480,6 +2480,125 @@ void RotationScaleMerge::RefineDecay(int n_groups) { // sit at the same s^2 but different frames, so the correction spreads them) - net harmful. slope = // 2 dB/dframe, so the total relative-B change across the run is |slope/2 * n_frames|. constexpr double DECAY_MIN_DELTA_B = 2.0; // A^2, minimum total relative-B over the run to engage + if (!segments.empty()) { + // Two hypotheses of how damage accrued over several sweeps, decided on held-out frames: a B rising + // linearly along the dose order (frame_dose), or one B per sweep - a step at each sweep boundary, + // flat within a sweep. Both are regressions of y = ln(I_ref / I_obs) on s^2 = 1/(4 d^2), the ramp + // on x = (dose - centre) s^2, the steps on x = s^2 with a slope per sweep, at the weights the slope + // fit above uses and within each frame (below). Each is fitted on the even frames and its + // weighted squared residual taken on the odd ones, against their own reference, and the other way + // round; the hypothesis that predicts the held-out frames better is applied, where it beats no + // decay at all by the margin a correction has to clear on one sweep (and the physical floor). + // The held-out frames' residual is read at weights fixed by the observations, so a correction + // cannot win by reshaping sigma. + using Sums = std::array; // Sw, Sx, Sy, Sxx, Sxy + struct Line { double a = 0.0, b = 0.0; }; + auto line_of = [](const Sums &S) { + Line l; + const double var = S[0] > 0.0 ? S[3] - S[1] * S[1] / S[0] : 0.0; + if (var > 0.0) l.b = (S[4] - S[1] * S[2] / S[0]) / var; + if (S[0] > 0.0) l.a = (S[2] - l.b * S[1]) / S[0]; + return l; + }; + auto reference = [&](int parity) { + std::fill(sw.begin(), sw.end(), 0.0); + std::fill(swI.begin(), swI.end(), 0.0); + for (const DecayTerm &o : t.term) { + if (parity >= 0 && (o.frame & 1) != parity) 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; + } + }; + // Every term of `parity` with its y and weight against the current reference. + auto for_terms = [&](int parity, const auto &fn) { + for (const DecayTerm &o : t.term) { + if ((parity >= 0 && (o.frame & 1) != parity) || 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; + fn(o, std::log(Iref / Is), (Is / sc) * (Is / sc)); + } + }; + auto x_ramp = [&](const DecayTerm &o) { return (dose(o) - fcenter) * s2_of(o.d); }; + // The terms of `parity` with y and both x taken about their frame's weighted mean: each frame has a + // scale of its own (the per-frame G), which took up its part of the B averaged over its + // reflections, so only how y varies with resolution WITHIN a frame is decay. + struct Point { int32_t seg; double w, y, xr, xs; }; + auto points = [&](int parity) { + reference(parity); + std::vector pts; + std::vector frame; + std::vector> mean(n_frames, std::array{}); // Sw, Sy, Sxr, Sxs + for_terms(parity, [&](const DecayTerm &o, double y, double w) { + pts.push_back({frame_segment[o.frame], w, y, x_ramp(o), s2_of(o.d)}); + frame.push_back(o.frame); + auto &m = mean[o.frame]; + m[0] += w; m[1] += w * y; m[2] += w * pts.back().xr; m[3] += w * pts.back().xs; + }); + for (size_t k = 0; k < pts.size(); ++k) { + const auto &m = mean[frame[k]]; + pts[k].y -= m[1] / m[0]; pts[k].xr -= m[2] / m[0]; pts[k].xs -= m[3] / m[0]; + } + return pts; + }; + struct Fit { double ramp = 0.0; std::vector step; }; + auto fit = [&](int parity) { + Sums r{}; + std::vector st(segments.size(), Sums{}); + for (const Point &q : points(parity)) { + r[0] += q.w; r[1] += q.w * q.xr; r[2] += q.w * q.y; r[3] += q.w * q.xr * q.xr; r[4] += q.w * q.xr * q.y; + Sums &a = st[q.seg]; + a[0] += q.w; a[1] += q.w * q.xs; a[2] += q.w * q.y; a[3] += q.w * q.xs * q.xs; a[4] += q.w * q.xs * q.y; + } + Fit out{line_of(r).b, {}}; + for (const Sums &a : st) out.step.push_back(line_of(a).b); + return out; + }; + // Held-out weighted squared residuals of no decay, the ramp and the steps. + auto score = [&](int test, const Fit &m, std::array &ssr) { + for (const Point &q : points(test)) { + const double e1 = q.y - m.ramp * q.xr, e2 = q.y - m.step[q.seg] * q.xs; + ssr[0] += q.w * q.y * q.y; ssr[1] += q.w * e1 * e1; ssr[2] += q.w * e2 * e2; + } + }; + std::array ssr{}; + const Fit even = fit(0), odd = fit(1); + score(1, even, ssr); + score(0, odd, ssr); + const Fit all = fit(-1); + + const double slope = all.ramp; // the slope decay_factor applies + std::vector b(segments.size()); + double b_mid = 0.0, b_lo = INFINITY, b_hi = -INFINITY; + for (size_t k = 0; k < segments.size(); ++k) { + b[k] = 0.5 * all.step[k]; + b_mid += b[k] * segments[k].n_frames / static_cast(n_frames); + b_lo = std::min(b_lo, b[k]); + b_hi = std::max(b_hi, b[k]); + } + const double delta_lin = std::fabs(0.5 * slope * n_frames), delta_step = b_hi - b_lo; + const bool lin_ok = delta_lin >= DECAY_MIN_DELTA_B && ssr[1] < (1.0 - CV_MIN_RELATIVE_GAIN) * ssr[0]; + const bool step_ok = delta_step >= DECAY_MIN_DELTA_B && ssr[2] < (1.0 - CV_MIN_RELATIVE_GAIN) * ssr[0]; + const bool use_step = step_ok && (!lin_ok || ssr[2] < ssr[1]); + const bool use_lin = lin_ok && !use_step; + segment_decay_model = use_step ? "per-sweep B" : use_lin ? "dose-linear" : "none"; + const double rel = ssr[0] > 0.0 ? 100.0 / ssr[0] : 0.0; + logger->Info("Decay over {} sweeps, held-out residual against no decay: dose-linear (dB = {:.2f} A^2 over " + "the dose order) {:.1f}%, one B per sweep (spread {:.2f} A^2) {:.1f}% => {}", segments.size(), + delta_lin, ssr[1] * rel, delta_step, ssr[2] * rel, segment_decay_model); + if (use_step || use_lin) + ParallelChunks(static_cast(fulls.size()), ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) { + Obs &o = fulls[i]; + if (!usable(o)) continue; + const double f = use_step ? std::exp(2.0 * (b[frame_segment[o.frame]] - b_mid) * s2_of(o.d)) + : decay_factor(o, slope); + o.corr = static_cast(o.corr * f); + } + }); + return; + } const double slope = fit_slope(-1); const double total_delta_B = std::fabs(0.5 * slope * n_frames); if (total_delta_B < DECAY_MIN_DELTA_B) { @@ -2505,25 +2624,22 @@ void RotationScaleMerge::RefineDecay(int n_groups) { 0.5 * slope * n_frames, 0.5 * slope); } -void RotationScaleMerge::MeasureSegmentDose(const DecayTerms &t, int n_groups) { - // A relative B per data set against the joint merge, as XSCALE fits it: Kabsch (2010) Acta Cryst. - // D66, 133-144; the dose read from it: Diederichs (2006) Acta Cryst. D62, 96-101. - // Radiation damage raises the relative B of a sweep collected later, so the B factors of the - // segments against their joint merge give the order the crystal took the dose in - read from the - // data, not from the order the sweeps were given in or the times their headers state. Per segment, - // the inverse-variance-of-log weighted slope of ln(I_ref / I_obs) against s^2 = 1/(4 d^2) is 2 B, - // with an intercept of its own (the scale is the per-frame G's). The least damaged segment is B = 0, - // and its frames start the dose axis. +std::vector RotationScaleMerge::SegmentB(const DecayTerms &t, int n_groups, int parity) const { + // Per segment, the inverse-variance-of-log weighted slope of ln(I_ref / I_obs) against + // s^2 = 1/(4 d^2), which is 2 B, with an intercept of its own (the scale is the per-frame G's), over + // the frames of `parity` (< 0: all) against a reference made of the same frames. NaN where a segment + // has nothing to fit. const size_t ns = segments.size(); std::vector sw(n_groups, 0.0), swI(n_groups, 0.0); for (const DecayTerm &o : t.term) { + if (parity >= 0 && (o.frame & 1) != parity) 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> S(ns, std::array{}); // Sw, Sx, Sy, Sxx, Sxy for (const DecayTerm &o : t.term) { - if (sw[o.group] <= 0.0) continue; + if ((parity >= 0 && (o.frame & 1) != parity) || 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; @@ -2539,6 +2655,20 @@ void RotationScaleMerge::MeasureSegmentDose(const DecayTerms &t, int n_groups) { if (var > 0.0) b[k] = 0.5 * (a[4] - a[1] * a[2] / a[0]) / var; } + return b; +} + +void RotationScaleMerge::MeasureSegmentDose(const DecayTerms &t, int n_groups) { + // A relative B per data set against the joint merge, as XSCALE fits it: Kabsch (2010) Acta Cryst. + // D66, 133-144; the dose read from it: Diederichs (2006) Acta Cryst. D62, 96-101. + // Radiation damage raises the relative B of a sweep collected later, so the B factors of the + // segments against their joint merge give the order the crystal took the dose in - read from the + // data, not from the order the sweeps were given in or the times their headers state. Per segment, + // the inverse-variance-of-log weighted slope of ln(I_ref / I_obs) against s^2 = 1/(4 d^2) is 2 B, + // with an intercept of its own (the scale is the per-frame G's). The least damaged segment is B = 0, + // and its frames start the dose axis. + const size_t ns = segments.size(); + const std::vector b = SegmentB(t, n_groups, -1); // A segment with nothing to measure keeps its place in the input order, after the measured ones. std::vector order(ns); std::iota(order.begin(), order.end(), 0); @@ -7489,6 +7619,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool refine_surfaces && (refine_decay_b || (along_frames && (absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0))); segment_results.assign(segments.size(), Result::Segment{}); + segment_decay_model.clear(); // 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. std::unique_ptr sub_mark; @@ -7602,6 +7733,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool segment_results[k].scale = median_of(g); } r.segments = segment_results; + r.segment_decay_model = segment_decay_model; // The fulls and the merge's working arrays are rebuilt by every arm before they are read, so they go // with it: kept, they sat beside the next arm's combine and through everything the caller did between // merges - on a large cell in P1, the fulls and about 330 bytes per ASU group, gigabytes. Faulting them diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 6b34a5129..f8886d385 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -157,6 +158,9 @@ public: int dose_rank = -1; // the segment's place in the dose order those B factors decide }; std::vector segments; + // How the decay over the segments was modelled - "dose-linear" (along the dose order), "per-sweep B" + // (one B a segment) or "none" - as the held-out frames decided; empty without segments or decay. + std::string segment_decay_model; }; // experiment: read live (its space group is changed by the caller between Run() calls). @@ -994,6 +998,9 @@ private: // Per segment: the relative B against the merged reference, and the dose order it puts the segments // in (segment_results, frame_dose). void MeasureSegmentDose(const DecayTerms &t, int n_groups); + // Each segment's B against a reference made of the same frames: those of `parity` (< 0: all). + std::vector SegmentB(const DecayTerms &t, int n_groups, int parity) const; + std::string segment_decay_model; // the last decay's choice (Result::segment_decay_model) // The terms the radiation-damage monitor built, handed on to RefineDecay: nothing between the two // touches corr, so they are the terms RefineDecay would build. Empty outside that stretch of RunArm. DecayTerms decay_terms_kept; diff --git a/tests/RotationScaleMergeSegmentsTest.cpp b/tests/RotationScaleMergeSegmentsTest.cpp index 1b33d0454..f20d05f49 100644 --- a/tests/RotationScaleMergeSegmentsTest.cpp +++ b/tests/RotationScaleMergeSegmentsTest.cpp @@ -21,7 +21,8 @@ namespace { struct SweepSpec { int n_frames; double scale; - double relative_b; + double relative_b; // at the sweep's first frame + double relative_b_end = NAN; // at its last; NaN: the same B throughout }; struct Synthetic { @@ -65,8 +66,11 @@ Synthetic MakeSweeps(const std::vector &sweeps, uint32_t seed) { for (int j = 0; j < 400; ++j) { const auto &m = hkl[pick(rng)]; const double s2 = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]) / (4.0 * CELL * CELL); + const double b = std::isfinite(sweep.relative_b_end) + ? sweep.relative_b + (sweep.relative_b_end - sweep.relative_b) * f / (sweep.n_frames - 1) + : sweep.relative_b; const double expected = out.I_true.at(FriedelKey(m[0], m[1], m[2])) * sweep.scale - * std::exp(-2.0 * sweep.relative_b * s2); + * std::exp(-2.0 * b * s2); const double sigma = std::sqrt(expected + 20.0); Reflection r{}; r.h = m[0]; r.k = m[1]; r.l = m[2]; @@ -170,6 +174,21 @@ TEST_CASE("RotationScaleMergeSegments_DoseOrderFromTheRelativeB", "[RotationScal CHECK(r.segments[1].relative_b == Catch::Approx(0.0).margin(1e-9)); CHECK(r.segments[0].relative_b == Catch::Approx(8.0).margin(1.0)); CHECK(r.segments[2].relative_b == Catch::Approx(4.0).margin(1.0)); + // Flat within each sweep: one B per sweep predicts the held-out frames better than a dose-linear ramp + CHECK(r.segment_decay_model == "per-sweep B"); +} + +TEST_CASE("RotationScaleMergeSegments_DoseLinearDecayAcrossSweeps", "[RotationScale]") { + // One crystal damaged steadily: the sweep given first was taken second, its B going on from where the + // other's ended. Along the dose order the B rises in a straight line, which the ramp describes. + const std::vector sweeps{{60, 1.0, 6.0, 12.0}, {60, 1.0, 0.0, 6.0}}; + auto data = MakeSweeps(sweeps, 15); + const auto r = Merge(data, SyntheticExperiment(true), Segments(sweeps)); + REQUIRE(r.segments.size() == 2); + CHECK(r.segments[1].dose_rank == 0); + CHECK(r.segments[0].dose_rank == 1); + CHECK(r.segments[0].relative_b == Catch::Approx(6.0).margin(1.5)); + CHECK(r.segment_decay_model == "dose-linear"); } TEST_CASE("RotationScaleMergeSegments_SegmentsMustCoverTheFrames", "[RotationScale]") {