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

RotationScaleMerge (tier F): decay / relative-B / radiation-damage sums over the fulls in fixed blocks

Browse Source 
Floating-point summation order only (tier F): the seven global sums of RefineDecay (slope fits,
held-out disagreements), the per-batch sums of FitRelativeBCurve, and the per-shell and
per-(batch, shell) sums of MeasureRadiationDamageB were single-threaded walks over every full.
They are now FixedBlockSum: the split depends on the number of fulls alone (ReductionBlocks,
ParallelBlocks), the block partials are merged in block order, so the sums have the same bits at
every thread count and on every machine - but not those of the serial walk.

Checked: p.mtz, p_P1.mtz, p.cif, p.hkl and the report identical to production on myob/cytc/thau
(GPU) and myob (CPU), and identical between -N 4/8/default (GPU cytc, myob; CPU myob). On these
sets the decay correction is not adopted (negligible or not cross-validated) and the monitor's
printed values do not move, so the F-tier change is not visible in their outputs; a set where the
decay correction is adopted will move in the last bits of corr.
Clean branch: p.mtz md5 identical on myob/cytc/thau/8a1a (GPU), cytc -N 8, myob -N 4 (CPU).
Measured (prototype with spans, 8a1a, 2 interleaved pairs, load 8-10): decay sections 2.70 -> 1.08 s,
radiation-damage monitor 1.04 -> 0.86 s, RSM total 37.1 -> 35.6 s.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
This commit is contained in:
leonarski_f
2026-10-07 06:30:20 +02:00
co-authored by Claude Opus 5.5
parent ea8132ceca
commit c8e2e6ef56
1 changed files with 66 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+66 -28
View File
@@ -239,6 +239,26 @@ namespace {
return out;
}
// A sum over [0, n) formed on all threads, in a split that depends on n alone: add(part, i) folds
// element i into its block's partial, and the partials are merged in block order (ParallelBlocks), so
// the sum has the same bits at every thread count - though not those of the one serial walk.
template <class S, class Add, class Merge>
S FixedBlockSum(int n, size_t nthreads, const S &zero, Add add, Merge merge) {
std::vector<S> part(ReductionBlocks(n), zero);
ParallelBlocks(n, nthreads, [&](int b, int lo, int hi) {
for (int i = lo; i < hi; ++i) add(part[b], i);
});
S total = zero;
for (const S &p : part) merge(total, p);
return total;
}
// Element-wise sum of equal-length arrays, the merge FixedBlockSum takes for per-batch / per-shell sums.
template <size_t K>
void AddArrays(std::array<std::vector<double>, K> &to, const std::array<std::vector<double>, K> &from) {
for (size_t k = 0; k < K; ++k)
for (size_t j = 0; j < to[k].size(); ++j) to[k][j] += from[k][j];
}
struct ScaleObs { double coeff, Iobs, weight; };
// Per-frame scale, linear in G: the weighted least-squares slope of the observations on their
@@ -2008,16 +2028,19 @@ void RotationScaleMerge::RefineDecay(int n_groups) {
sw[g] = s_w; swI[g] = s_wI;
}
});
double Sw = 0, Sx = 0, Sy = 0, Sxx = 0, Sxy = 0;
for (const DecayTerm &o : t.term) {
if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) continue;
using Lin = std::array<double, 5>; // Sw, Sx, Sy, Sxx, Sxy
const Lin S = FixedBlockSum(static_cast<int>(t.term.size()), nthreads, Lin{},
[&](Lin &a, int k) {
const DecayTerm &o = t.term[k];
if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) return;
const double Iref = swI[o.group] / sw[o.group];
const double Is = static_cast<double>(o.I) * o.corr, sc = static_cast<double>(o.sigma) * o.corr;
if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue;
if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) return;
const double w = (Is / sc) * (Is / sc);
const double x = (o.image_number - fcenter) * s2_of(o.d), y = std::log(Iref / Is);
Sw += w; Sx += w * x; Sy += w * y; Sxx += w * x * x; Sxy += w * x * y;
}
a[0] += w; a[1] += w * x; a[2] += w * y; a[3] += w * x * x; a[4] += w * x * y;
}, [](Lin &a, const Lin &b) { for (int j = 0; j < 5; ++j) a[j] += b[j]; });
const double Sw = S[0], Sx = S[1], Sy = S[2], Sxx = S[3], Sxy = S[4];
const double var = (Sw > 0.0) ? Sxx - Sx * Sx / Sw : 0.0;
// Clamp only to guard against a near-collinear (var ~ 0) blow-up; generous enough to reach very
// strong damage (slope = 2 dB/dframe, so +-1.0 admits total relative-B up to ~n_frames/2 A^2).
@@ -2045,14 +2068,16 @@ void RotationScaleMerge::RefineDecay(int n_groups) {
sw[g] = s_w; swI[g] = s_wI;
}
});
double num = 0.0, den = 0.0;
for (const DecayTerm &o : t.term) {
if ((o.frame & 1) != parity || sw[o.group] <= 0.0) continue;
using ND = std::array<double, 2>;
const ND nd = FixedBlockSum(static_cast<int>(t.term.size()), nthreads, ND{}, [&](ND &a, int k) {
const DecayTerm &o = t.term[k];
if ((o.frame & 1) != parity || sw[o.group] <= 0.0) return;
const double Iref = swI[o.group] / sw[o.group];
const double Is = static_cast<double>(o.I) * o.corr * decay_factor(o, slope);
if (!std::isfinite(Iref) || Iref <= 0.0) continue;
num += std::abs(Is - Iref); den += Iref;
}
if (!std::isfinite(Iref) || Iref <= 0.0) return;
a[0] += std::abs(Is - Iref); a[1] += Iref;
}, [](ND &a, const ND &b) { a[0] += b[0]; a[1] += b[1]; });
const double num = nd[0], den = nd[1];
return den > 0.0 ? num / den : 0.0;
};
@@ -2158,16 +2183,20 @@ std::vector<double> RotationScaleMerge::FitRelativeBCurve(const DecayTerms &t, i
sw[g] = s_w; swI[g] = s_wI;
}
});
std::vector<double> num(n_batch, 0.0), den(n_batch, 0.0);
for (const DecayTerm &o : t.term) {
if ((gparity >= 0 && (o.group & 1) != gparity) || sw[o.group] <= 0.0) continue;
using Batches = std::array<std::vector<double>, 2>; // num, den per batch
const Batches nd = FixedBlockSum(static_cast<int>(t.term.size()), nthreads,
Batches{std::vector<double>(n_batch, 0.0), std::vector<double>(n_batch, 0.0)},
[&](Batches &a, int k) {
const DecayTerm &o = t.term[k];
if ((gparity >= 0 && (o.group & 1) != gparity) || sw[o.group] <= 0.0) return;
const double Iref = swI[o.group] / sw[o.group];
const double Is = static_cast<double>(o.I) * o.corr, sc = static_cast<double>(o.sigma) * o.corr;
if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue;
if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) return;
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;
}
a[0][batch_of(o)] += w * s2 * y; a[1][batch_of(o)] += w * s2 * s2;
}, AddArrays<2>);
const std::vector<double> &num = nd[0], &den = nd[1];
std::vector<double> b = SolveCurvatureSmoothedB(num, den, RELATIVE_B_MAX);
double dw = 0.0, dbw = 0.0;
for (int c = 0; c < n_batch; ++c) { dw += den[c]; dbw += den[c] * b[c]; }
@@ -2253,12 +2282,16 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) {
// spread over all of it would spend most of the grid on noise, leaving a batch with too few to fit at
// all. Pooled over the whole run a shell either has signal or it has none, so walk out to the first
// that has none and lay the shells inside that.
std::vector<double> rw(MONITOR_SHELLS, 0.0), rwI(MONITOR_SHELLS, 0.0);
for (const DecayTerm &o : t.term) {
using Shells = std::array<std::vector<double>, 2>; // rw, rwI per shell
const Shells r_shell = FixedBlockSum(static_cast<int>(t.term.size()), nthreads,
Shells{std::vector<double>(MONITOR_SHELLS, 0.0), std::vector<double>(MONITOR_SHELLS, 0.0)},
[&](Shells &a, int j) {
const DecayTerm &o = t.term[j];
const double sc = static_cast<double>(o.sigma) * o.corr, w = 1.0 / (sc * sc);
const int k = shell_of(s2_of(o.d));
rw[k] += w; rwI[k] += w * static_cast<double>(o.I) * o.corr;
}
a[0][k] += w; a[1][k] += w * static_cast<double>(o.I) * o.corr;
}, AddArrays<2>);
const std::vector<double> &rw = r_shell[0], &rwI = r_shell[1];
int n_live = 0;
while (n_live < MONITOR_SHELLS && rw[n_live] > 0.0
&& rwI[n_live] / std::sqrt(rw[n_live]) >= MONITOR_MIN_SHELL_ISIGMA)
@@ -2278,15 +2311,20 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) {
// mean is well determined where a single observation is not, admits negative intensities, and carries
// its own I/sigma - which is what says whether the batch can be measured at all.
const int n_cell = n_batch * MONITOR_SHELLS;
std::vector<double> sw(n_cell, 0.0), swI(n_cell, 0.0), swR(n_cell, 0.0), sws2(n_cell, 0.0);
for (const DecayTerm &o : t.term) {
if (sw0[o.group] <= 0.0) continue;
using Cells = std::array<std::vector<double>, 4>; // sw, swI, swR, sws2 per (batch, shell)
const std::vector<double> zero_cells(n_cell, 0.0);
const Cells cells = FixedBlockSum(static_cast<int>(t.term.size()), nthreads,
Cells{zero_cells, zero_cells, zero_cells, zero_cells},
[&](Cells &a, int j) {
const DecayTerm &o = t.term[j];
if (sw0[o.group] <= 0.0) return;
const double Iref = swI0[o.group] / sw0[o.group];
if (!std::isfinite(Iref)) continue;
if (!std::isfinite(Iref)) return;
const double s2 = s2_of(o.d), sc = static_cast<double>(o.sigma) * o.corr, w = 1.0 / (sc * sc);
const int i = batch_of(o) * MONITOR_SHELLS + shell_of(s2);
sw[i] += w; swI[i] += w * static_cast<double>(o.I) * o.corr; swR[i] += w * Iref; sws2[i] += w * s2;
}
a[0][i] += w; a[1][i] += w * static_cast<double>(o.I) * o.corr; a[2][i] += w * Iref; a[3][i] += w * s2;
}, AddArrays<4>);
const std::vector<double> &sw = cells[0], &swI = cells[1], &swR = cells[2], &sws2 = cells[3];
// Per batch: slope of ln(<I_ref0> / <I_obs>) against s^2 over the shells that hold signal, weighted by
// the shell's own (I/sigma)^2 since that is the precision of its log. Fitted WITH an intercept and the
// normal equations centred, so that only the resolution-DEPENDENT part reaches the curve: a batch that