mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 90 Wiki Activity

Rotation scaling: frame guards and delta-CC1/2 read the frames by what they carry

Browse Source 
Three follow-ups to the converging scale loop, all on sweeps with a stretch the crystal barely
diffracted on:

- The sweep ledger's scale channel (and with it delta-CC1/2's "normal frame" test) and the
  space-group search's scale floor measure a frame against the precision-weighted typical frame
  G_ref (TypicalFrameScale) instead of the run median: on a sweep that spent most of its turn out
  of the beam the median frame is itself a dead one, and every stretch then reads as typical. The
  two collapse guards stay on the median on purpose - a frame that collapsed toward zero and was
  not quite dropped has its fulls re-fitted with a scale of 1/G, a cubic mean is then theirs, and
  the floor read against it dropped three quarters of every live frame's fulls (measured).

- delta-CC1/2's sigma-tau statistic enters each reflection with the information it carries (the
  same per-frame factor as the CC1/2 weight), so a dead stretch, whose scaled-up noise flooded the
  mean error variance of every mixed reflection and read as harm up to the 25% cap, now costs
  about nothing and is left to the ledger.

- The scale loop's step test weighs each frame by its merge weight (its observations at its scale
  squared): a dead frame's scale is fitted on noise and wanders by orders of magnitude every
  iteration, carries nothing into the merge, is dropped after the loop, and must not hold the loop
  open.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
This commit is contained in:
leonarski_fandClaude Fable 5.1 committed 2026-09-20 18:45:17 +02:00
1 parent 004839a8ee
commit ddc7855900
1 file changed
+66 -27
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+66 -27
View File
@@ -1430,15 +1430,23 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop(
}
rescale(scaled_any, ratio);
if (it > 0) {
double ss = 0.0;
// The step, weighted by what each frame carries into the merge - its observations at its
// scale squared, which is how the inverse-variance weight of a scaled observation goes.
// A frame the crystal did not diffract on has its scale fitted on noise and wanders by
// orders of magnitude every iteration (measured on a sweep 60% out of the beam: scales
// at 1e-7 of the typical frame, still moving at the cap); it carries nothing into the
// merge and is dropped from it after the loop, so it does not hold the loop open either.
double ss = 0.0, sw = 0.0;
int n = 0;
for (int f = 0; f < n_frames; ++f)
if (scaled_any[f] && g_prev[f] > 0.0 && g[f] > 0.0) {
const double d = std::log(g[f] / g_prev[f]);
ss += d * d;
const double w = frame_obs_count[f] * g[f] * g[f];
ss += w * d * d;
sw += w;
++n;
}
out.step = n > 0 ? std::sqrt(ss / n) : 0.0;
out.step = sw > 0.0 ? std::sqrt(ss / sw) : 0.0;
logger.Debug("Per-frame scaling ({}): iteration {}, rms |dlogG| {:.2e} over {} frames, "
"gauge {:.4f}", what, it + 1, out.step, n, g_ref);
if (out.step < SCALING_TOLERANCE) {
@@ -1959,21 +1967,30 @@ void RotationScaleMerge::MeasureSweepQuality(const std::vector<uint8_t> &partial
return;
// Run gauges. Only the ratio to the run matters - the absolute scale is degenerate with the merge's.
std::vector<double> g_fitted, cc_fitted;
// The scale's gauge is the precision-weighted typical frame (TypicalFrameScale), not the median
// frame: on a sweep that spent most of its turn out of the beam the median frame IS a dead one,
// and every stretch then reads as typical (measured: a run whose good frames sat at 1.3-1.7 and
// whose median was 0.21 reported its dead stretch at "scale 0.24"). The CC channel keeps its
// median - a correlation is not a scale, and dead frames carry no CC to drag it with.
std::vector<double> cc_fitted;
std::vector<int64_t> n_obs(n_frames, 0);
int n_scaled = 0;
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 (partial_scaled[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0) {
n_obs[f] = frame_count[f];
++n_scaled;
}
if (cc_n[f] >= MIN_REFLECTIONS_FOR_IMAGE_CC && std::isfinite(cc[f]))
cc_fitted.push_back(cc[f]);
}
if (static_cast<int>(g_fitted.size()) < SWEEP_MIN_SCALED_FRAMES)
if (n_scaled < SWEEP_MIN_SCALED_FRAMES)
return;
auto median_of = [](std::vector<double> &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 g_typ = TypicalFrameScale(g_partial, n_obs);
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
@@ -2177,15 +2194,24 @@ namespace {
// the variance of the unique reflections' mean intensities against the error variance of a half
// dataset's mean, accumulated one reflection at a time. No half-set split, so no random numbers -
// the same answer every run, and more precise than the split estimate.
//
// Each reflection enters with the information it carries: the precision its observations would
// have at the run's typical frame scale over the precision they have (the same per-frame factor
// as the CC1/2 weight, CCHalfFrameFactors), so a reflection measured only on frames the crystal
// barely diffracted on, whose scaled-up noise would otherwise flood the mean error variance of
// every mixed reflection, counts for what it knows. Unweighted, a dead stretch read as harm
// (-0.15 on a sweep 60% out of the beam) and was convicted up to the cap; weighted, it costs
// about nothing and is left to the ledger, which is what it is - a stretch that adds nothing.
struct SigmaTau {
double n = 0, sum_mean = 0, sum_mean2 = 0, sum_eps = 0;
void Add(double mean, double eps) {
n += 1; sum_mean += mean; sum_mean2 += mean * mean; sum_eps += eps;
double n = 0, w = 0, sum_mean = 0, sum_mean2 = 0, sum_eps = 0;
void Add(double mean, double eps, double weight) {
n += 1; w += weight; sum_mean += weight * mean; sum_mean2 += weight * mean * mean;
sum_eps += weight * eps;
}
[[nodiscard]] double CC() const {
if (n < 4) return NAN;
const double var_y = (sum_mean2 - sum_mean * sum_mean / n) / (n - 1);
const double eps = sum_eps / n;
if (n < 4 || !(w > 0.0)) return NAN;
const double var_y = (sum_mean2 - sum_mean * sum_mean / w) / w * n / (n - 1);
const double eps = sum_eps / w;
return var_y + eps / 2 > 0.0 ? (var_y - eps / 2) / (var_y + eps / 2) : NAN;
}
};
@@ -2293,12 +2319,17 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) {
// The observations, once: scaled intensity, ASU group (-1 = not in this merge, the same test the
// merge itself applies), and the frame order that makes any frame range a contiguous span.
std::vector<double> obs_I(n_full);
std::vector<double> obs_I(n_full), obs_info(n_full);
std::vector<int32_t> obs_group(n_full), order(n_full), fstart(n_frames + 1, 0);
std::vector<int64_t> n_obs_frame(n_frames, 0);
for (int i = 0; i < n_full; ++i)
if (UsableFull(fulls[i])) ++n_obs_frame[fulls[i].frame];
const std::vector<double> frame_factor = CCHalfFrameFactors(frame_scale_fulls, n_obs_frame);
for (int i = 0; i < n_full; ++i) {
const Obs &o = fulls[i];
obs_I[i] = o.I * o.corr;
obs_group[i] = UsableFull(o) ? o.group : -1;
obs_info[i] = 1.0 / frame_factor[o.frame]; // the information the observation carries, <= 1
++fstart[o.frame + 1];
}
for (int f = 0; f < n_frames; ++f) fstart[f + 1] += fstart[f];
@@ -2310,16 +2341,18 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) {
// Per-reflection totals over the frames still in, and the per-reflection contribution of one frame
// range - which is all a leave-one-out needs, because a count, a sum and a sum of squares subtract.
std::vector<int32_t> tot_n(n_groups, 0), d_n(n_groups, 0), touched;
std::vector<double> tot_s1(n_groups, 0.0), tot_s2(n_groups, 0.0),
d_s1(n_groups, 0.0), d_s2(n_groups, 0.0);
std::vector<double> tot_s1(n_groups, 0.0), tot_s2(n_groups, 0.0), tot_w(n_groups, 0.0),
d_s1(n_groups, 0.0), d_s2(n_groups, 0.0), d_w(n_groups, 0.0);
auto build_totals = [&]() {
std::fill(tot_n.begin(), tot_n.end(), 0);
std::fill(tot_s1.begin(), tot_s1.end(), 0.0);
std::fill(tot_s2.begin(), tot_s2.end(), 0.0);
std::fill(tot_w.begin(), tot_w.end(), 0.0);
for (int i = 0; i < n_full; ++i) {
const int g = obs_group[i];
if (g < 0 || frame_rejected[fulls[i].frame]) continue;
tot_n[g] += 1; tot_s1[g] += obs_I[i]; tot_s2[g] += obs_I[i] * obs_I[i];
tot_w[g] += obs_info[i];
}
};
@@ -2340,6 +2373,7 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) {
if (g < 0) continue;
if (d_n[g] == 0) touched.push_back(g);
d_n[g] += 1; d_s1[g] += obs_I[i]; d_s2[g] += obs_I[i] * obs_I[i];
d_w[g] += obs_info[i];
}
}
// Only reflections that are still measurable WITHOUT the range count, which is eq. 5's own
@@ -2351,11 +2385,11 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) {
if (tot_n[g] >= 2 && n_without >= 2) {
double mean, eps;
ReflectionMoments(tot_n[g], tot_s1[g], tot_s2[g], mean, eps);
with.Add(mean, eps);
with.Add(mean, eps, tot_w[g] / tot_n[g]);
ReflectionMoments(n_without, tot_s1[g] - d_s1[g], tot_s2[g] - d_s2[g], mean, eps);
without.Add(mean, eps);
without.Add(mean, eps, (tot_w[g] - d_w[g]) / n_without);
}
d_n[g] = 0; d_s1[g] = 0.0; d_s2[g] = 0.0;
d_n[g] = 0; d_s1[g] = 0.0; d_s2[g] = 0.0; d_w[g] = 0.0;
}
const double cc_with = with.CC(), cc_without = without.CC();
if (with.n < DELTA_CC_MIN_REFLECTIONS || !std::isfinite(cc_with) || !std::isfinite(cc_without))
@@ -3180,6 +3214,11 @@ bool RotationScaleMerge::DropCollapsedScales(const std::vector<uint8_t> &fitted_
apply.assign(n, 0);
ratio.assign(n, 1.0);
// Against the MEDIAN frame, not the precision-weighted typical frame the sweep ledger reads: a
// frame that collapsed toward zero and was not quite dropped has its intensities amplified by
// 1/G, the fulls built from them are re-fitted with a scale of that size, and a cubic mean is then
// theirs - measured on a sweep 60% out of the beam, the fulls' precision-weighted typical scale
// sat on such frames and the floor dropped three quarters of every live frame's fulls.
std::vector<double> fitted;
fitted.reserve(n);
for (int f = 0; f < n; ++f)
@@ -3256,6 +3295,7 @@ bool RotationScaleMerge::DropCollapsedFullScales(bool from_staging) {
&& o.corr != 1.0f)
g_frame[o.frame] = 1.0 / static_cast<double>(o.corr);
// The median, not the precision-weighted typical frame - see DropCollapsedScales.
std::vector<double> fitted;
fitted.reserve(g_frame.size());
for (const double gf : g_frame)
@@ -4834,14 +4874,13 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
// --- 2a'. Search pass only: drop the frames the crystal barely diffracted on, whose 1/G would
// amplify their intensities into a P1 merge with no multiplicity to catch them. ---
if (for_search) {
std::vector<double> fitted;
std::vector<int64_t> n_obs(n_frames, 0);
for (int f = 0; f < n_frames; ++f)
if (frame_scaled_scratch[f] && std::isfinite(g_partial[f]) && g_partial[f] > 0.0)
fitted.push_back(g_partial[f]);
if (!fitted.empty()) {
const size_t mid = fitted.size() / 2;
std::nth_element(fitted.begin(), fitted.begin() + mid, fitted.end());
const double floor_g = fitted[mid] * SEARCH_MIN_SCALE_RATIO;
n_obs[f] = frame_count[f];
if (std::any_of(n_obs.begin(), n_obs.end(), [](int64_t c) { return c > 0; })) {
// Against the precision-weighted typical frame - see MeasureSweepQuality.
const double floor_g = TypicalFrameScale(g_partial, n_obs) * SEARCH_MIN_SCALE_RATIO;
std::vector<uint8_t> reject(n_frames, 0);
int n_rejected = 0;
for (int f = 0; f < n_frames; ++f)
@@ -4865,7 +4904,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
if (reject[partials[i].frame]) partials[i].corr = 0.0f;
});
logger.Info("Space-group search: ignoring {} of {} frames whose scale came out below "
"1/{:.0f} of the run median (the crystal barely diffracted there)",
"1/{:.0f} of the run's typical frame (the crystal barely diffracted there)",
n_rejected, n_frames, 1.0 / SEARCH_MIN_SCALE_RATIO);
}
}