diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 6c5b11e6..fe58bd68 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -25,6 +25,7 @@ #include "../../common/CrystalLattice.h" #include "../../common/Definitions.h" #include "../../common/JFJochException.h" +#include "../../common/JFJochMath.h" #include "../../common/ResolutionShells.h" namespace { @@ -44,6 +45,29 @@ namespace { // took the merged CC1/2 from 93.7% to 17.2%). 0.02 is 3.5x below anything real and 12x above the // failure. constexpr double MIN_CREDIBLE_SCALE_RATIO = 0.02; + + // --- Sweep-quality diagnostic (MeasureSweepQuality) --- + // A stretch is reported only when BOTH per-frame channels are down: the scale (how much the crystal + // diffracted) and the CC to merge (whether what it diffracted is still usable). The CC channel is what + // keeps a merely attenuated stretch out - absorption and flux scale a frame's intensities without + // changing how well they correlate with the merged reference, and a strong crystal seen through a 4x + // absorption dip is still good data. + constexpr int SWEEP_MIN_SCALED_FRAMES = 20; // fewer than this and the run gauges mean nothing + constexpr double SWEEP_WINDOW_DEG = 5.0; // running-median window for both channels + constexpr double SWEEP_MIN_RANGE_DEG = 5.0; // shorter than this is a fluctuation, not a stretch + constexpr double SWEEP_SCALE_FRACTION = 0.5; // scale below this fraction of the run median + constexpr double SWEEP_CC_FRACTION = 0.7; // ... and CC below this fraction of the run median + constexpr double SWEEP_DEAD_FRACTION = 0.25; // scaled frames below this: nothing was recorded + constexpr double SWEEP_LOST_FRACTION = 0.8; // scaled frames below this x the run's: frames lost + constexpr double SWEEP_DECAY_CC_FRACTION = 0.7; // CC already this far down before a terminal range + constexpr double SWEEP_FULL_TURN_DEG = 350.0; // a once-per-revolution claim needs a whole revolution + constexpr double SWEEP_HARMONIC_R2 = 0.7; // the fundamental must explain this much of log-scale + constexpr double SWEEP_HARMONIC_RATIO = 1.5; // ... and dominate the 180 deg (crystal-shape) term + constexpr double SWEEP_HARMONIC_DEPTH = 1.5; // ... at this peak-to-trough of the FITTED fundamental + // (the observed curve swings further; the gate that + // guarantees the modulation costs data is the CC dip) + constexpr double SWEEP_HARMONIC_CC_DIP = 0.8; // ... and cost real signal at its trough + constexpr float MAX_FRAME_GAP = 2.0f; // a rocking event is a run of frames no more apart than this constexpr double CHI2_1_MEDIAN = 0.454936; // A post-scale-fulls correction surface (decay / absorption) is applied only if its held-out @@ -1134,6 +1158,244 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { for (int c = 0; c < n_batch; ++c) rad_damage_b_batch[c] = static_cast(b[c]); } +namespace { + // Running median over `window` frames. A median rather than a mean because the real per-frame scale + // moves fast (its median frame-to-frame step in log is ~0.06) and one dead frame inside a good stretch + // must not drag the window down. + std::vector RunningMedian(const std::vector &v, int window) { + const int n = static_cast(v.size()), half = window / 2; + std::vector out(n), buf; + for (int i = 0; i < n; ++i) { + buf.assign(v.begin() + std::max(0, i - half), v.begin() + std::min(n, i + half + 1)); + const size_t mid = buf.size() / 2; + std::nth_element(buf.begin(), buf.begin() + mid, buf.end()); + out[i] = buf[mid]; + } + return out; + } + + // Least-squares fit of one cos(k*t)/sin(k*t) pair to `y` over the frames flagged in `use`, subtracted + // from y in place and added into `model`. Returns the amplitude. Over a full turn the harmonics are + // near-orthogonal, so fitting them one after another gives the joint answer without a matrix solve. + double RemoveHarmonic(std::vector &y, std::vector &model, + const std::vector &use, int k) { + const int n = static_cast(y.size()); + double scc = 0, sss = 0, scs = 0, syc = 0, sys = 0; + for (int f = 0; f < n; ++f) { + if (!use[f]) continue; + const double t = 2.0 * PI * k * f / n, c = std::cos(t), s = std::sin(t); + scc += c * c; sss += s * s; scs += c * s; syc += y[f] * c; sys += y[f] * s; + } + const double det = scc * sss - scs * scs; + if (!(std::fabs(det) > 1e-12)) + return 0.0; + const double a = (syc * sss - sys * scs) / det, b = (sys * scc - syc * scs) / det; + for (int f = 0; f < n; ++f) { + const double t = 2.0 * PI * k * f / n, fit = a * std::cos(t) + b * std::sin(t); + y[f] -= fit; + model[f] += fit; + } + return std::hypot(a, b); + } +} + +void RotationScaleMerge::MeasureSweepQuality(const std::vector &partial_scaled, + const std::vector &cc, + const std::vector &cc_n) { + sweep_quality = SweepQuality{}; + const auto gon = x.GetGoniometer(); + const double osc_deg = gon ? std::fabs(gon->GetIncrement_deg()) : 0.0; + if (n_frames < 2 || !(osc_deg > 1e-6)) + return; + + // Run gauges. Only the ratio to the run matters - the absolute scale is degenerate with the merge's. + std::vector g_fitted, cc_fitted; + for (int f = 0; f < n_frames; ++f) { + if (partial_scaled[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0) + g_fitted.push_back(g_partial[f]); + if (cc_n[f] >= MIN_REFLECTIONS_FOR_IMAGE_CC && std::isfinite(cc[f])) + cc_fitted.push_back(cc[f]); + } + if (static_cast(g_fitted.size()) < SWEEP_MIN_SCALED_FRAMES) + return; + auto median_of = [](std::vector &v) { + const size_t mid = v.size() / 2; + std::nth_element(v.begin(), v.begin() + mid, v.end()); + return v[mid]; + }; + const double g_typ = median_of(g_fitted); + const double cc_typ = cc_fitted.empty() ? 0.0 : median_of(cc_fitted); + + // The two channels. The scale is the sample side alone: DivideOutIncidentFlux took the beam out of the + // observations before G was ever fitted, so what is left in G is the crystal - volume in the beam, + // absorption, or damage. A frame that got no scale contributed nothing and reads 0 in both channels. + // Without a per-frame CC the protection against reporting a merely attenuated stretch is gone, so + // leave that channel at 1 and report nothing rather than report it on the scale alone. + const bool have_cc = cc_typ > 0.0; + std::vector scale(n_frames, 0.0), quality(n_frames, have_cc ? 0.0 : 1.0); + for (int f = 0; f < n_frames; ++f) { + if (partial_scaled[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0) + scale[f] = g_partial[f] / g_typ; + if (have_cc && cc_n[f] >= MIN_REFLECTIONS_FOR_IMAGE_CC && std::isfinite(cc[f])) + quality[f] = cc[f] / cc_typ; + } + + sweep_quality.measured = true; + sweep_quality.sweep_deg = static_cast(n_frames * osc_deg); + { // How much the beam itself moved, for scale: this is the part already divided out of the above. + // 5th to 95th percentile, not min to max - a single frame whose background was measured off a + // handful of reflections would otherwise set a headline number. + std::vector flux; + for (int f = 0; f < n_frames; ++f) + if (frame_flux[f] > 0.0) flux.push_back(frame_flux[f]); + if (flux.size() >= 20) { + std::sort(flux.begin(), flux.end()); + const double lo = flux[flux.size() / 20], hi = flux[flux.size() - 1 - flux.size() / 20]; + if (lo > 0.0) sweep_quality.flux_peak_to_trough = static_cast(hi / lo); + } + } + + // One cycle of modulation per revolution. A crystal off the rotation axis leaves the illuminated + // volume once per turn; the crystal's own shape absorbs on a 180 deg period, so the fundamental is + // what separates the two - Evans (Acta Cryst. D62 (2006) 72-82) notes that illuminated volume and + // absorption are otherwise indistinguishable. Fitted on the log scale after a linear trend has taken + // the dose out, and only over a sweep long enough to have seen the crystal come back. + // + // This one test runs on the TOTAL scale, flux and all, unlike everything else here. The flux proxy is + // the mean background, and a crystal drifting out of the illuminated volume takes its own diffuse and + // solvent scattering with it, so dividing that out removes part of the very modulation being looked + // for. The beam, on the other hand, cannot be periodic in the goniometer angle - it does not know + // where the goniometer is - so a once-per-turn component of the background belongs to the sample. + bool modulated = false; + int trough = 0, mod_first = 0, mod_last = -1; + if (sweep_quality.sweep_deg >= SWEEP_FULL_TURN_DEG) { + std::vector use(n_frames, 0); + std::vector y(n_frames, 0.0), model(n_frames, 0.0); + double sw = 0, sf = 0, sy = 0, sff = 0, sfy = 0; + for (int f = 0; f < n_frames; ++f) { + if (!(scale[f] > 0.0) || !(frame_flux[f] > 0.0)) continue; + use[f] = 1; + y[f] = std::log(scale[f] * frame_flux[f]); + sw += 1; sf += f; sy += y[f]; sff += double(f) * f; sfy += double(f) * y[f]; + } + const double det = sw * sff - sf * sf; + const double slope = det > 0.0 ? (sw * sfy - sf * sy) / det : 0.0; + const double icept = sw > 0.0 ? (sy - slope * sf) / sw : 0.0; + double ss_tot = 0.0, ss_lin = 0.0; + const double mean_y = sw > 0.0 ? sy / sw : 0.0; + for (int f = 0; f < n_frames; ++f) { + model[f] = icept + slope * f; + if (!use[f]) continue; + ss_tot += (y[f] - mean_y) * (y[f] - mean_y); + y[f] -= model[f]; + ss_lin += y[f] * y[f]; + } + const double amp1 = RemoveHarmonic(y, model, use, 1); + double ss_h1 = 0.0; + for (int f = 0; f < n_frames; ++f) if (use[f]) ss_h1 += y[f] * y[f]; + const double amp2 = RemoveHarmonic(y, model, use, 2); + const double r2_gain = ss_tot > 0.0 ? (ss_lin - ss_h1) / ss_tot : 0.0; + const double depth = std::exp(2.0 * amp1); + // model now holds trend + both harmonics; the fundamental's trough is where the crystal is worst. + int peak = 0; + for (int f = 1; f < n_frames; ++f) { + if (model[f] < model[trough]) trough = f; + if (model[f] > model[peak]) peak = f; + } + const auto qm = RunningMedian(quality, std::max(3, static_cast(std::lround(SWEEP_WINDOW_DEG / osc_deg)))); + // A modulation is a dimming, not a disappearance: the crystal has to be measurably worse at the + // trough than at the peak, and still delivering data there. A sweep whose trough is simply dead + // is a sweep with a dead arc, and is reported as one. + modulated = r2_gain >= SWEEP_HARMONIC_R2 && amp1 >= SWEEP_HARMONIC_RATIO * amp2 + && depth >= SWEEP_HARMONIC_DEPTH && qm[peak] > 0.0 && qm[trough] > 0.0 + && qm[trough] < SWEEP_HARMONIC_CC_DIP * qm[peak]; + if (modulated) { + sweep_quality.modulation_peak_to_trough = static_cast(depth); + // Report the bottom quarter of the modulation around its trough: the part of the turn where + // the loss is worst, not the whole half-cycle that is merely below average. + const double lim = model[trough] + 0.25 * (model[peak] - model[trough]); + mod_first = mod_last = trough; + while (mod_first > 0 && model[mod_first - 1] < lim) --mod_first; + while (mod_last + 1 < n_frames && model[mod_last + 1] < lim) ++mod_last; + sweep_quality.ranges.push_back({mod_first, mod_last, SweepQualityReason::LossOfCentring}); + } + } + + // Contiguous stretches where BOTH channels are down. Gaps shorter than the minimum range are closed: + // a stretch interrupted by a few good frames is one event, not two. + const int window = std::max(3, static_cast(std::lround(SWEEP_WINDOW_DEG / osc_deg))); + const int min_len = std::max(1, static_cast(std::lround(SWEEP_MIN_RANGE_DEG / osc_deg))); + const auto sm = RunningMedian(scale, window); + const auto qm = RunningMedian(quality, window); + std::vector> found; + for (int f = 0; f < n_frames; ++f) { + if (!(sm[f] < SWEEP_SCALE_FRACTION && qm[f] < SWEEP_CC_FRACTION)) continue; + if (!found.empty() && f - found.back().second - 1 <= min_len) found.back().second = f; + else found.emplace_back(f, f); + } + + for (const auto &[first, last] : found) { + if (last - first + 1 < min_len) + continue; + // The modulation trough, if there is one, is already reported; do not report it twice. + if (first <= mod_last && mod_first <= last) + continue; + sweep_quality.ranges.push_back({first, last, SweepQualityReason::WeakDiffraction}); + } + std::sort(sweep_quality.ranges.begin(), sweep_quality.ranges.end(), + [](const SweepQualityRange &a, const SweepQualityRange &b) { return a.first_image < b.first_image; }); + + // Fill in each range's numbers, and say what it is. The order is from the most specific evidence to + // the least: nothing recorded at all; then a decay that had already set in before the range and runs + // to the end of the sweep; then the once-per-revolution modulation; then whether frames were lost or + // only intensity. What is left is a loss of diffracting power whose cause these data do not fix. + double run_scaled = 0.0; + for (int f = 0; f < n_frames; ++f) run_scaled += scale[f] > 0.0 ? 1.0 : 0.0; + run_scaled /= n_frames; + + for (auto &r : sweep_quality.ranges) { + const int n_in = r.last_image - r.first_image + 1; + double sum_s = 0, sum_q = 0, n_ok = 0, sum_b = 0, n_b = 0; + for (int f = r.first_image; f <= r.last_image; ++f) { + sum_s += scale[f]; + sum_q += quality[f]; + n_ok += scale[f] > 0.0 ? 1.0 : 0.0; + if (!rad_damage_b_batch.empty() && rad_damage_batch_deg > 0.0) { + const int c = std::min(rad_damage_b_batch.size() - 1, + static_cast(f * osc_deg / rad_damage_batch_deg)); + sum_b += rad_damage_b_batch[c]; n_b += 1; + } + } + r.rotation_deg = static_cast(n_in * osc_deg); + r.mean_relative_scale = static_cast(sum_s / n_in); + r.mean_relative_cc = static_cast(sum_q / n_in); + r.indexed_fraction = static_cast(n_ok / n_in); + r.severity = static_cast(std::clamp(1.0 - sum_s / n_in, 0.0, 1.0)); + r.relative_b = n_b > 0 ? static_cast(sum_b / n_b) : NAN; + if (r.reason == SweepQualityReason::LossOfCentring) + continue; + // Radiation damage is progressive: the per-frame CC has to have been falling BEFORE the range, + // and the range has to run to the end of the sweep. A crystal that simply leaves the beam at the + // end fails the first test, and a crystal that recovers fails the second. + bool decayed = false; + if (r.last_image >= n_frames - std::max(1, min_len / 2) && r.first_image > n_frames / 5) { + const int k = std::max(1, r.first_image / 10); + double first_tenth = 0, last_tenth = 0; + for (int f = 0; f < k; ++f) first_tenth += quality[f]; + for (int f = r.first_image - k; f < r.first_image; ++f) last_tenth += quality[f]; + decayed = first_tenth > 0.0 && last_tenth < SWEEP_DECAY_CC_FRACTION * first_tenth; + } + if (r.indexed_fraction < SWEEP_DEAD_FRACTION) + r.reason = SweepQualityReason::NoDiffraction; + else if (decayed) + r.reason = SweepQualityReason::RadiationDamage; + else if (r.indexed_fraction < SWEEP_LOST_FRACTION * run_scaled) + r.reason = SweepQualityReason::CrystalOutOfBeam; + else + r.reason = SweepQualityReason::WeakDiffraction; + } +} + 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 @@ -2426,6 +2688,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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; + out.sweep_quality = sweep_quality; // 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). @@ -2724,6 +2987,10 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { // (report-only; captures the full damage signature, not a residual). Skipped on the P1 search pass. if (!for_search) MeasureRadiationDamageB(n_groups); + // Sweep-quality diagnostic, on the per-frame scale the partial scaling just fitted (the flux is + // already out of it) and the per-frame CC computed above. Report-only; drops nothing. + if (!for_search) + MeasureSweepQuality(partial_scaled, cc, cc_n); if (!for_search && refine_decay_b) RefineDecay(n_groups); if (!for_search && relative_b_deg > 0.0) diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 8e8ea824..9e0345b8 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -185,6 +185,10 @@ private: 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) + // Sweep-quality diagnostic (MeasureSweepQuality; report-only, copied into the result statistics by + // MergeAndStats). Empty and not measured until it runs. + SweepQuality sweep_quality; + // Working per-group arrays (sized to the current group count; reused). std::vector group_h, group_k, group_l; @@ -269,6 +273,11 @@ private: // 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); + // Sweep-quality diagnostic (report-only): find the contiguous stretches of the sweep over which the + // crystal delivered much less than the rest of the run, and say what each one looks like. Reads the + // per-frame scale (with the incident flux already divided out) and the per-frame CC to merge. + void MeasureSweepQuality(const std::vector &partial_scaled, const std::vector &cc, + const std::vector &cc_n); // 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 diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 0c7599cf..05a3dcf2 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2072,6 +2072,34 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("{}", os.str()); result.radiation_damage_text = os.str(); } + + // Sweep-quality report: the stretches of the sweep over which the crystal delivered much + // less than the rest of the run, and what each one looks like. Nothing is excluded because + // of it - the frames still carry signal, and this is a message for the beamline, not a + // filter. Frame numbers are processed-image ordinals, inclusive, as in _image.dat. + const auto &sq = sm.statistics.sweep_quality; + if (sq.measured) { + std::ostringstream os; + os << fmt::format("Sweep quality over {:.0f} deg (per-image scale with the incident flux, " + "which varied {:.2f}x, already divided out):\n", + sq.sweep_deg, sq.flux_peak_to_trough); + if (sq.ranges.empty()) { + os << " no stretch of the sweep is materially worse than the run"; + } else { + os << " frames rotation diagnosis severity scale CC indexed\n"; + for (const auto &r : sq.ranges) + os << fmt::format(" {:<17s} {:6.1f} deg {:<16s} {:5.2f} {:5.2f} {:5.2f} {:4.0f}%\n", + fmt::format("{}-{}", r.first_image, r.last_image), r.rotation_deg, + SweepQualityReasonText(r.reason), r.severity, r.mean_relative_scale, + r.mean_relative_cc, 100.0 * r.indexed_fraction); + os << " => severity is the fraction of the run's typical diffracting power missing " + "over the range"; + if (sq.modulation_peak_to_trough >= 1.05f) + os << fmt::format("\n => once-per-revolution modulation of the per-image scale: " + "{:.1f}x peak to trough", sq.modulation_peak_to_trough); + } + logger.Info("{}", os.str()); + } } if (result.consensus_cell && write_files && config_.write_merged) {