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]") {