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

Merge cpu-scaling-tail into rc175: CPU scaling and merging walks on all threads in the same order, fewer bytes moved per iteration, output unchanged and independent of the thread count

Browse Source 
--no-gpu Scaling and merging (provisional, loaded machine): 8qaw 120 -> 77 s, 9v16 71 -> 41 s, 4nwv 74 -> 57 s.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-15 22:48:01 +02:00
2 parents 496145a399 4a8aff41c2
commit 0a8eda37fe
2 files changed
+401 -204
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/scale_merge/RotationScaleMerge.cpp
+334 -185
View File
@@ -1009,6 +1009,7 @@ void RotationScaleMerge::BuildInRangeObservations(std::unique_ptr<SortKey[]> &ke
// the product that is measured from the data rather than computed from the geometry.
obs.prescaling_corr = static_cast<float>(r.prescaling_corr * r.qe_corr * r.flight_corr / frame_flux[o]);
obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg;
obs.var_per_I = 0.0f; obs.capture = 0.0f;
obs.px = r.predicted_x; obs.py = r.predicted_y;
obs.image_number = r.image_number;
obs.frame = o;
@@ -1461,11 +1462,12 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) {
return n_groups;
}
#endif
// Without the GPU: stamp the ASU-group id per obs from its raw hkl + the precomputed finiteness, on
// the partials and in a flat array the group CSR below is counted from.
// Without the GPU: stamp the ASU-group id per obs from its raw hkl + the precomputed finiteness, in
// the flat array the CPU passes read it from (partial_group) and the group CSR below is counted from.
// Plain arrays, here and for gperm below: every element is written, and a sized vector would
// first zero all of them on one thread - once per merge, a sweep's worth of observations.
std::unique_ptr<int32_t[]> group_ids(new int32_t[n_partials_obs]);
partial_group.reset(new int32_t[n_partials_obs]);
int32_t *group_ids = partial_group.get();
// Parallel over raw-hkl runs: distinct runs own disjoint perm ranges, hence disjoint observations.
ParallelChunks(n_run, nthreads, [&](int rlo, int rhi) {
for (int r = rlo; r < rhi; ++r) {
@@ -1473,9 +1475,7 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) {
const int lo = rawrun_start[r], hi = rawrun_start[r] + rawrun_count[r];
for (int p = lo; p < hi; ++p) {
const int i = perm[p];
const int gi = (g >= 0 && finite_ok[i]) ? g : -1;
group_ids[i] = gi;
partials[i].group = gi;
group_ids[i] = (g >= 0 && finite_ok[i]) ? g : -1;
}
}
});
@@ -1562,7 +1562,8 @@ void RotationScaleMerge::ReducePartialGroupMeans(int n_groups, std::vector<doubl
const int lo = partial_group_start[g], hi = lo + partial_group_count[g];
for (int p = lo; p < hi; ++p) {
float I_corr = 0.0f;
const double w = ReferenceWeight(partials[partial_group_perm[p]], min_partiality, I_corr);
const Obs &o = partials[partial_group_perm[p]];
const double w = ReferenceWeight(g, o.corr, o.partiality, o.I, o.sigma, min_partiality, I_corr);
if (w <= 0.0) continue;
sw += w;
swI += w * I_corr;
@@ -1573,9 +1574,9 @@ void RotationScaleMerge::ReducePartialGroupMeans(int n_groups, std::vector<doubl
}
void RotationScaleMerge::ReduceScalingGroupMeans(int n_groups, const int32_t *group_start,
const int32_t *group_count, const std::vector<ScalingRef> &refs,
const std::vector<float> &corr,
std::vector<double> &out_mean) const {
const int32_t *group_count, const FillVector<ScalingRef> &refs,
FillVector<float> &corr,
std::vector<double> &out_mean, const CorrUpdate *update) const {
out_mean.resize(n_groups);
ParallelChunks(n_groups, ThreadsForWork(corr.size(), nthreads), [&](int glo, int ghi) {
for (int g = glo; g < ghi; ++g) {
@@ -1583,6 +1584,7 @@ void RotationScaleMerge::ReduceScalingGroupMeans(int n_groups, const int32_t *gr
const int lo = group_start[g], hi = lo + group_count[g];
for (int p = lo; p < hi; ++p) {
const ScalingRef &r = refs[p];
if (update) update->Apply(r, corr[p]);
float I_corr = 0.0f;
const double w = ReferenceWeight(g, corr[p], r.partiality, r.I, r.sigma, min_partiality, I_corr);
if (w <= 0.0) continue;
@@ -1669,16 +1671,21 @@ double RotationScaleMerge::ChooseLogScaleSmoothing(const std::vector<double> &y,
return best_lambda;
}
void RotationScaleMerge::FitScalingFrames(const std::vector<ScalingFit> &fit, const std::vector<int32_t> &fit_start,
void RotationScaleMerge::FitScalingFrames(const FillVector<ScalingFit> &fit, const std::vector<int32_t> &fit_start,
size_t min_obs, const std::vector<double> &group_mean_in,
std::vector<double> &g) {
const int n_fr = static_cast<int>(fit_start.size()) - 1;
const int n_fit = static_cast<int>(fit.size());
constexpr int FIT_PREFETCH = 16;
std::vector<uint8_t> scaled(n_fr, 0);
std::vector<double> info(n_fr, 0.0);
ParallelFor(n_fr, nthreads, [&](int f) {
double num = 0.0, den = 0.0;
size_t n = 0;
for (int i = fit_start[f]; i < fit_start[f + 1]; ++i) {
// The group mean is a gather from an array as long as the group count - the cost of this
// loop on a P1 merge - so it is asked for FIT_PREFETCH observations ahead.
if (i + FIT_PREFETCH < n_fit) Prefetch(&group_mean_in[fit[i + FIT_PREFETCH].group]);
const ScalingFit &o = fit[i];
// The reference is the group's mean WITH this observation in it. Without it the update is
// an undamped Jacobi step on the sparse graph a P1 merge of a short sweep is, and it
@@ -1689,9 +1696,11 @@ void RotationScaleMerge::FitScalingFrames(const std::vector<ScalingFit> &fit, co
// iterations, which is what RunScalingLoop's tolerance is for.
const double mean = group_mean_in[o.group];
if (!std::isfinite(mean)) continue;
const double coeff = o.coeff_per_mean * mean;
num += o.w2 * coeff * o.Iobs;
den += o.w2 * coeff * coeff;
const double weight = SafeInv(o.sigma, 1.0), w2 = weight * weight;
const double coeff_per_mean = o.partiality * SafeInv(o.prescaling_corr, 1.0);
const double coeff = coeff_per_mean * mean;
num += w2 * coeff * static_cast<double>(o.I);
den += w2 * coeff * coeff;
++n;
}
if (n < min_obs) return;
@@ -3740,6 +3749,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
// Terms per block of the per-cell sums in the fit below: every block zeroes and hands back two
// ncell-long arrays, so a block has to hold many more terms than there are cells.
constexpr int SURFACE_BLOCK = 32768;
constexpr int SURFACE_PREFETCH = 16;
bool on_gpu = false;
int n_terms = 0;
#ifdef JFJOCH_USE_CUDA
@@ -3889,24 +3899,33 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
// terms in fulls order, because the counting sort is stable (it walks `term`, which is in fulls
// order). A separate copy rather than a permutation: chasing every term through an index at random
// over a few hundred megabytes is the locality this pass used to stay serial to keep.
struct RefTerm { float I, sigma, corr; int32_t cell; uint8_t parity; };
std::vector<int32_t> g_start;
std::vector<RefTerm> gterm;
// One copy per frame parity: most reference passes are of one half (the cross-validation), and
// each then reads only its own. A pass over all the terms walks a group's two halves together in
// term order (`at`, the term's position), which is the fulls order the single copy held.
struct RefTerm { float I, sigma, corr; int32_t cell, at; };
std::array<std::vector<int32_t>, 2> g_start;
std::array<std::vector<RefTerm>, 2> gterm;
if (!on_gpu) {
g_start.assign(n_groups + 1, 0);
for (const Term &t : term) ++g_start[t.group + 1];
for (int g = 0; g < n_groups; ++g) g_start[g + 1] += g_start[g];
gterm.resize(term.size());
std::vector<int32_t> fill(g_start.begin(), g_start.end() - 1);
for (int par = 0; par < 2; ++par) g_start[par].assign(n_groups + 1, 0);
for (size_t k = 0; k < term.size(); ++k) ++g_start[term_parity[k]][term[k].group + 1];
std::array<std::vector<int32_t>, 2> fill;
for (int par = 0; par < 2; ++par) {
for (int g = 0; g < n_groups; ++g) g_start[par][g + 1] += g_start[par][g];
gterm[par].resize(g_start[par][n_groups]);
fill[par].assign(g_start[par].begin(), g_start[par].end() - 1);
}
for (size_t k = 0; k < term.size(); ++k) {
const Term &t = term[k];
gterm[fill[t.group]++] = {t.I, t.sigma, t.corr, t.cell, term_parity[k]};
const int par = term_parity[k];
gterm[par][fill[par][t.group]++] = {t.I, t.sigma, t.corr, t.cell, static_cast<int32_t>(k)};
}
}
// The per-group reference of one subset with surface A applied - the first pass of both the fit
// and the score.
std::vector<double> sw(n_groups), swI(n_groups);
// iref is the reference itself, swI / sw, NaN where a group has no weight: the fit below reads one
// array at random per term instead of two (CPU only).
std::vector<double> sw(n_groups), swI(n_groups), iref(on_gpu ? 0 : n_groups);
auto reference = [&](int parity, const std::vector<double> &A) {
#ifdef JFJOCH_USE_CUDA
if (on_gpu) {
@@ -3917,9 +3936,7 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
ParallelChunks(n_groups, nt, [&](int glo, int ghi) {
for (int g = glo; g < ghi; ++g) {
double s_w = 0.0, s_wI = 0.0;
for (int k = g_start[g]; k < g_start[g + 1]; ++k) {
const RefTerm &t = gterm[k];
if (parity >= 0 && t.parity != parity) continue;
auto add = [&](const RefTerm &t) {
// Keep the intensity's own product together before it is weighted: w * ((I*corr)*a)
// is what the loops this replaced summed, and w * I * corr * a rounds differently.
// The score decides whether a surface is kept, so a last-place difference here can
@@ -3927,8 +3944,17 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
const double a = A[t.cell], Is = static_cast<double>(t.I) * t.corr * a;
const double sc = static_cast<double>(t.sigma) * t.corr * a, w = 1.0 / (sc * sc);
s_w += w; s_wI += w * Is;
};
if (parity >= 0) {
for (int k = g_start[parity][g]; k < g_start[parity][g + 1]; ++k) add(gterm[parity][k]);
} else {
int k0 = g_start[0][g], k1 = g_start[1][g];
const int e0 = g_start[0][g + 1], e1 = g_start[1][g + 1];
while (k0 < e0 || k1 < e1)
add(k1 == e1 || (k0 < e0 && gterm[0][k0].at < gterm[1][k1].at) ? gterm[0][k0++] : gterm[1][k1++]);
}
sw[g] = s_w; swI[g] = s_wI;
iref[g] = s_w > 0.0 ? s_wI / s_w : std::numeric_limits<double>::quiet_NaN();
}
});
};
@@ -3991,9 +4017,11 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector<int32_t> &cell, int
std::fill(xcross.begin(), xcross.end(), 0.0);
std::fill(xref2.begin(), xref2.end(), 0.0);
for (int k = lo; k < hi; ++k) {
// The reference is a gather from an array as long as the group count, so it is
// asked for SURFACE_PREFETCH terms ahead.
if (k + SURFACE_PREFETCH < hi) Prefetch(&iref[term[sel[k + SURFACE_PREFETCH]].group]);
const Term &o = term[sel[k]];
if (sw[o.group] <= 0.0) continue;
const double Iref = swI[o.group] / sw[o.group], a = A[o.cell];
const double Iref = iref[o.group], a = A[o.cell]; // NaN: the group has no weight
const double Is = static_cast<double>(o.I) * o.corr * a, sc = static_cast<double>(o.sigma) * o.corr * a;
if (!std::isfinite(Iref) || Iref <= 0.0 || !(sc > 0.0)) continue;
const double w = 1.0 / (sc * sc);
@@ -4352,7 +4380,8 @@ bool RotationScaleMerge::UsableFull(const Obs &o) const {
}
void RotationScaleMerge::Combine() {
fulls.clear();
// `fulls` is not cleared: SortFullsByFrame writes every record, so its old contents only save the
// resize there from constructing them again.
g_full.assign(n_frames, 1.0);
fulls_dropped_overloaded = 0;
const int n_run = static_cast<int>(rawrun_start.size());
@@ -4514,6 +4543,7 @@ void RotationScaleMerge::Combine() {
};
size_t n_used = 0;
std::vector<std::vector<Obs>> part; // the fulls of consecutive runs, in run order
if (!observation_dump_path.empty() || nthreads <= 1) {
// Serial (diagnostic dump needs a single writer).
@@ -4524,15 +4554,16 @@ void RotationScaleMerge::Combine() {
// is a per-full batch (ZD) and px/py the predicted detector position (XD/YD).
dump << "# h k l I sigma d n_frames captured_fraction peak_frame px py\n";
}
fulls.reserve(n_run);
part.resize(1);
part[0].reserve(n_run);
for (int r = 0; r < n_run; ++r)
n_used += process_rawrun(r, fulls, dump.is_open() ? &dump : nullptr, fulls_dropped_overloaded);
n_used += process_rawrun(r, part[0], dump.is_open() ? &dump : nullptr, fulls_dropped_overloaded);
} else {
// Parallel over contiguous rawrun chunks; concatenate the per-thread fulls in run order so the
// result is deterministic.
// Parallel over contiguous rawrun chunks; the per-thread fulls are taken in run order below, so
// the result is deterministic.
const int nt = static_cast<int>(std::min(nthreads, static_cast<size_t>(n_run)));
const int chunk = (n_run + nt - 1) / nt;
std::vector<std::vector<Obs>> part(nt);
part.resize(nt);
std::vector<size_t> used(nt, 0);
std::vector<int64_t> overloaded(nt, 0);
std::vector<std::future<void>> futures;
@@ -4547,44 +4578,74 @@ void RotationScaleMerge::Combine() {
}));
}
for (auto &f : futures) f.get();
size_t total = 0;
for (int t = 0; t < nt; ++t) {
total += part[t].size(); n_used += used[t]; fulls_dropped_overloaded += overloaded[t];
n_used += used[t]; fulls_dropped_overloaded += overloaded[t];
}
fulls.reserve(total);
for (int t = 0; t < nt; ++t)
fulls.insert(fulls.end(), part[t].begin(), part[t].end());
}
SortFullsByFrame();
SortFullsByFrame(part);
logger->Info("3D combine: {} fulls from {} partials, {} rocking events dropped for a saturated pixel",
fulls.size(), n_used, fulls_dropped_overloaded);
}
void RotationScaleMerge::BuildFullsGroupCSR(int n_groups) {
// Each group's fulls in index order. (group, index) is a total order, so the sort on all threads
// returns the one sequence the serial counting sort filled.
const int nf = static_cast<int>(fulls.size());
fulls_group_start.assign(n_groups + 1, 0);
for (int i = 0; i < nf; ++i)
if (fulls[i].group >= 0) ++fulls_group_start[fulls[i].group + 1];
for (int g = 0; g < n_groups; ++g) fulls_group_start[g + 1] += fulls_group_start[g];
fulls_group_perm.resize(fulls_group_start[n_groups]);
std::vector<int32_t> fill(fulls_group_start.begin(), fulls_group_start.end() - 1);
for (int i = 0; i < nf; ++i)
if (fulls[i].group >= 0) fulls_group_perm[fill[fulls[i].group]++] = i;
const std::vector<int32_t> grouped = KeptIndices(nf, nthreads, [&](int i) { return fulls[i].group >= 0; });
const int n = static_cast<int>(grouped.size());
std::vector<std::pair<int32_t, int32_t>> key(n);
ParallelChunks(n, ThreadsForWork(key.size(), nthreads), [&](int lo, int hi) {
for (int j = lo; j < hi; ++j) key[j] = {fulls[grouped[j]].group, grouped[j]};
});
ParallelSort(key.begin(), key.end(), nthreads, std::less<std::pair<int32_t, int32_t>>());
fulls_group_start.resize(n_groups + 1);
fulls_group_perm.resize(n);
// Group g starts where the first key of a group >= g sits, so a group's start is written by the key
// that opens it, together with the empty groups just below it.
ParallelChunks(n, ThreadsForWork(key.size(), nthreads), [&](int lo, int hi) {
for (int j = lo; j < hi; ++j) {
fulls_group_perm[j] = key[j].second;
const int g = key[j].first;
if (j > 0 && key[j - 1].first == g) continue;
for (int e = j > 0 ? key[j - 1].first + 1 : 0; e <= g; ++e) fulls_group_start[e] = j;
}
});
for (int e = n > 0 ? key[n - 1].first + 1 : 0; e <= n_groups; ++e) fulls_group_start[e] = n;
}
void RotationScaleMerge::SortFullsByFrame() {
// Sort the fulls by their (peak) frame and build per-frame CSR ranges, so the scale-fulls step can
// fit a per-frame G by slicing contiguous ranges (the same layout the partials use).
std::sort(fulls.begin(), fulls.end(),
[](const Obs &a, const Obs &b) { return a.frame < b.frame; });
void RotationScaleMerge::SortFullsByFrame(const std::vector<std::vector<Obs>> &part) {
// Sort the fulls - the parts one after another - by their (peak) frame and build per-frame CSR
// ranges, so the scale-fulls step can fit a per-frame G by slicing contiguous ranges (the same layout
// the partials use). What std::sort does to a sequence depends on nothing but the outcomes of the
// comparisons it makes, so sorting each full's (frame, position) under the same frame-only test
// leaves the positions in exactly the order the eighty-byte records would have been left in - ties
// included - while moving eight bytes instead of eighty. The records are then gathered once.
struct FrameAt { int32_t frame, at; };
std::vector<size_t> part_start(part.size() + 1, 0);
for (size_t t = 0; t < part.size(); ++t) part_start[t + 1] = part_start[t] + part[t].size();
const int nf = static_cast<int>(part_start.back());
std::vector<FrameAt> order(nf);
ParallelFor(static_cast<int>(part.size()), nthreads, [&](int t) {
for (size_t j = 0; j < part[t].size(); ++j)
order[part_start[t] + j] = {part[t][j].frame, static_cast<int32_t>(part_start[t] + j)};
});
std::sort(order.begin(), order.end(), [](const FrameAt &a, const FrameAt &b) { return a.frame < b.frame; });
fulls.resize(nf);
ParallelChunks(nf, ThreadsForWork(order.size(), nthreads), [&](int lo, int hi) {
for (int i = lo; i < hi; ++i) {
const size_t at = order[i].at;
const size_t t = std::upper_bound(part_start.begin(), part_start.end(), at) - part_start.begin() - 1;
fulls[i] = part[t][at - part_start[t]];
}
});
fulls_frame_start.assign(n_frames, 0);
fulls_frame_count.assign(n_frames, 0);
for (int i = 0; i < static_cast<int>(fulls.size()); ) {
const int f = fulls[i].frame;
for (int i = 0; i < nf; ) {
const int f = order[i].frame;
fulls_frame_start[f] = i;
int j = i;
while (j < static_cast<int>(fulls.size()) && fulls[j].frame == f) ++j;
while (j < nf && order[j].frame == f) ++j;
fulls_frame_count[f] = j - i;
i = j;
}
@@ -4601,12 +4662,13 @@ void RotationScaleMerge::ComputePerFrameCC(const std::vector<double> &partial_gr
const int lo = frame_start[f], hi = frame_start[f] + frame_count[f];
for (int i = lo; i < hi; ++i) {
const auto &o = partials[i];
const int g = partial_group[i];
if (o.on_ice) continue;
if (o.group < 0) continue;
if (g < 0) continue;
if (o.partiality < min_partiality) continue;
if (!std::isfinite(o.I) || !std::isfinite(o.corr) || o.corr <= 0.0f) continue;
if (!std::isfinite(o.sigma) || o.sigma <= 0.0f) continue;
const double mean = partial_group_mean[o.group];
const double mean = partial_group_mean[g];
if (!std::isfinite(mean)) continue;
const double img = static_cast<double>(o.I) * o.corr;
sx += img; sy += mean; sx2 += img * img; sy2 += mean * mean; sxy += img * mean;
@@ -4852,23 +4914,34 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
for (int g = 0; g < n_groups; ++g)
if (cnt[g] >= 2 && sw[g] > 0.0) em_mean[g] = swI[g] / sw[g];
samples.clear();
samples.reserve(fulls.size());
for (const auto &o : fulls) {
if (!usable_merge(o) || cnt[o.group] < 2) continue;
const size_t s = n_hand ? 2 * static_cast<size_t>(o.group) + obs_hand[&o - fulls.data()] : 0;
const bool on_hand = n_hand && group_has_hands[o.group] && cnth[s] >= 2 && swh[s] > 0.0;
const double mean = on_hand ? swIh[s] / swh[s] : em_mean[o.group];
if (!std::isfinite(mean)) continue;
// One sample per full that gives one, in index order: which fulls those are, then each one's
// sample on its own - on all threads.
auto sample_of = [&](int i, Sample &out) {
const int g = mf.group[i]; // usable_merge
if (g < 0 || cnt[g] < 2) return false;
const Obs &o = fulls[i];
const size_t s = n_hand ? 2 * static_cast<size_t>(g) + obs_hand[i] : 0;
const bool on_hand = n_hand && group_has_hands[g] && cnth[s] >= 2 && swh[s] > 0.0;
const double mean = on_hand ? swIh[s] / swh[s] : em_mean[g];
if (!std::isfinite(mean)) return false;
const double sigma_corr = static_cast<double>(o.sigma) * o.corr;
const double s2 = sigma_corr * sigma_corr;
const double w = 1.0 / s2;
const double factor = 1.0 - w / (on_hand ? swh[s] : sw[o.group]);
if (factor < 0.05) continue;
const double factor = 1.0 - w / (on_hand ? swh[s] : sw[g]);
if (factor < 0.05) return false;
const double resid = static_cast<double>(o.I) * o.corr - mean;
samples.push_back({counting_variance(o, em_mean[o.group], s2), mean * mean,
resid * resid / factor, o.d});
}
out = {counting_variance(o, em_mean[g], s2), mean * mean, resid * resid / factor, o.d};
return true;
};
const std::vector<int32_t> sampled = KeptIndices(n_full, nthreads, [&](int i) {
Sample unused;
return sample_of(i, unused);
});
samples.resize(sampled.size());
ParallelChunks(static_cast<int>(sampled.size()), ThreadsForWork(sampled.size(), nthreads),
[&](int lo, int hi) {
for (int j = lo; j < hi; ++j) sample_of(sampled[j], samples[j]);
});
}
// Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on
@@ -5339,42 +5412,53 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
did_gpu_acc = true;
}
#endif
if (!did_gpu_acc)
for (const auto &o : fulls) {
if (!usable_merge(o)) continue;
if (rejected_obs[&o - fulls.data()]) { // by the Wilson test above
++reject_count;
continue;
}
const float I_corr = o.I * o.corr;
float sigma_corr = o.sigma * o.corr;
sigma_corr = corrected_sigma(o, I_corr, sigma_corr);
if (reject_outliers && error_model_active && std::isfinite(reject_median[o.group])) {
// The pooled pair, widened by the shell's measured anomalous variance.
const double v_add = reject_var_add.empty() ? 0.0 : reject_var_add[o.group];
const double I_b = std::isfinite(em_mean[o.group]) ? em_mean[o.group] : I_corr;
if (OutsideOutlierBand(I_corr, reject_median[o.group], reject_nsigma,
static_cast<double>(sigma_corr) * sigma_corr + v_add, I_b,
error_model_b * I_b, reject_up, reject_down)) {
++reject_count;
rejected_obs[&o - fulls.data()] = 1;
continue;
// The threads split the GROUPS, each walking a group's fulls over the group CSR, so every group
// is still summed in index order; the rejection count is an integer sum, added up per block.
if (!did_gpu_acc) {
std::vector<size_t> blk_rejected(ReductionBlocks(n_groups, 512), 0);
ParallelBlocks(n_groups, nthreads, [&](int blk, int g0, int g1) {
for (int g = g0; g < g1; ++g) {
for (int p = fulls_group_start[g]; p < fulls_group_start[g + 1]; ++p) {
const int i = fulls_group_perm[p];
if (mf.group[i] != g) continue; // usable_merge
const Obs &o = fulls[i];
if (rejected_obs[i]) { // by the Wilson test above
++blk_rejected[blk];
continue;
}
const float I_corr = o.I * o.corr;
float sigma_corr = o.sigma * o.corr;
sigma_corr = corrected_sigma(o, I_corr, sigma_corr);
if (reject_outliers && error_model_active && std::isfinite(reject_median[o.group])) {
// The pooled pair, widened by the shell's measured anomalous variance.
const double v_add = reject_var_add.empty() ? 0.0 : reject_var_add[o.group];
const double I_b = std::isfinite(em_mean[o.group]) ? em_mean[o.group] : I_corr;
if (OutsideOutlierBand(I_corr, reject_median[o.group], reject_nsigma,
static_cast<double>(sigma_corr) * sigma_corr + v_add, I_b,
error_model_b * I_b, reject_up, reject_down)) {
++blk_rejected[blk];
rejected_obs[i] = 1;
continue;
}
}
const double w = 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
const double wI = w * I_corr;
auto &a = acc[o.group];
const int half = merge_half[i];
a.swI += wI; a.sw += w;
a.swIh[half] += wI; a.swh[half] += w; a.nh[half]++;
a.swh_typ[half] += w * frame_cc_factor[o.frame];
// Any contributing observation on an ice ring marks the merged reflection (as
// MergeOnTheFly does on the stills path): the intensity is still merged and written,
// the flag only lets a consumer that must not read ring signal as crystal signal -
// the CC1/2 resolution-cutoff fit - leave it out.
if (o.on_ice) a.on_ice = true;
if (!std::isfinite(a.d) && std::isfinite(o.d) && o.d > 0.0f) a.d = o.d;
}
}
const double w = 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
const double wI = w * I_corr;
auto &a = acc[o.group];
const int half = merge_half[&o - fulls.data()];
a.swI += wI; a.sw += w;
a.swIh[half] += wI; a.swh[half] += w; a.nh[half]++;
a.swh_typ[half] += w * frame_cc_factor[o.frame];
// Any contributing observation on an ice ring marks the merged reflection (as
// MergeOnTheFly does on the stills path): the intensity is still merged and written,
// the flag only lets a consumer that must not read ring signal as crystal signal -
// the CC1/2 resolution-cutoff fit - leave it out.
if (o.on_ice) a.on_ice = true;
if (!std::isfinite(a.d) && std::isfinite(o.d) && o.d > 0.0f) a.d = o.d;
}
}, 512);
for (size_t r : blk_rejected) reject_count += r;
}
// A merged reflection's d from the cell the completeness denominator is enumerated with, not
// from whichever frame first observed it: a per-frame cell puts a reflection near a shell edge
// on the other side of it from its possible twin, and a shell then reads above 100 %.
@@ -5765,31 +5849,45 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
did_gpu_rmeas = true;
}
#endif
if (!did_gpu_rmeas)
for (const auto &o : fulls) {
if (o.group < 0 || acc[o.group].d < d_floor) continue;
if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge -> also out of R_meas
if (!frame_cell_ok[o.frame]) continue;
if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue;
if (o.partiality < min_partiality) continue;
const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr;
if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue;
const auto shell = shells.GetShell(o.d);
if (!shell || *shell < 0 || *shell >= n_shells) continue;
// Only what the merge kept counts towards multiplicity - see the GPU branch above. The
// R_meas accumulation below keeps its own, wider filter.
if (usable_merge(o))
sa[*shell].total_obs++;
if (std::isfinite(merged_I[o.group])) {
auto &r = rmeas[o.group];
const double sc = corrected_sigma(o, I_corr, sigma_corr);
const double v = 1.0 / (sc * sc);
const double dev = std::fabs(static_cast<double>(I_corr) - merged_I[o.group]);
r.sum_abs_dev += dev; r.sum_I += I_corr;
r.w_abs_dev += v * dev; r.w_I += v * I_corr; r.sum_v += v; r.sum_v2 += v * v;
r.n++; r.shell = *shell;
if (!did_gpu_rmeas) {
// The threads split the GROUPS, each walking a group's fulls over the group CSR, so every
// group is still summed in index order; the observation counts are integers, so the
// blocks may add theirs up in any split.
const int nb = ReductionBlocks(n_groups, 512);
std::vector<std::vector<int>> blk_obs(nb, std::vector<int>(n_shells, 0));
ParallelBlocks(n_groups, nthreads, [&](int blk, int g0, int g1) {
for (int g = g0; g < g1; ++g) {
if (acc[g].d < d_floor) continue;
for (int p = fulls_group_start[g]; p < fulls_group_start[g + 1]; ++p) {
const int i = fulls_group_perm[p];
const Obs &o = fulls[i];
if (rejected_obs[i]) continue; // outlier-rejected in the merge -> also out of R_meas
if (!frame_cell_ok[o.frame]) continue;
if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue;
if (o.partiality < min_partiality) continue;
const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr;
if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue;
const auto shell = shells.GetShell(o.d);
if (!shell || *shell < 0 || *shell >= n_shells) continue;
// Only what the merge kept counts towards multiplicity - see the GPU branch above.
// The R_meas accumulation below keeps its own, wider filter.
if (mf.group[i] == g) // usable_merge
blk_obs[blk][*shell]++;
if (std::isfinite(merged_I[g])) {
auto &r = rmeas[g];
const double sc = corrected_sigma(o, I_corr, sigma_corr);
const double v = 1.0 / (sc * sc);
const double dev = std::fabs(static_cast<double>(I_corr) - merged_I[g]);
r.sum_abs_dev += dev; r.sum_I += I_corr;
r.w_abs_dev += v * dev; r.w_I += v * I_corr; r.sum_v += v; r.sum_v2 += v * v;
r.n++; r.shell = *shell;
}
}
}
}
}, 512);
for (const auto &c : blk_obs)
for (int sh = 0; sh < n_shells; ++sh) sa[sh].total_obs += c[sh];
}
for (const auto &r : rmeas) {
if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue;
double num, den;
@@ -5871,6 +5969,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
pair_hkl[pair_of_group[g]] = group_pair[g];
}
std::vector<AnomAcc> anom(pair_hkl.size());
if (own_table) anom_export.reserve(anom.size());
// The half-sets CCanom needs are balanced within the MATE, not merely within the reflection:
// it correlates one half's I(+)-I(-) against the other's, so each hand separately needs both
// of its halves filled, and a mate left entirely in one half loses the whole pair rather than
@@ -5949,11 +6048,19 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
std::vector<float> row_mean(anom.size() * 2, NAN);
std::vector<int32_t> row_shell(anom.size() * 2, -1), row_of_cell(anom.size() * 2);
for (size_t c = 0; c < row_of_cell.size(); ++c) row_of_cell[c] = static_cast<int32_t>(c);
// Whether each pair is centric, on its own, so on all threads: the test walks the group's
// operations, and it is the costly part of the serial pair walk below.
std::vector<uint8_t> pair_centric(anom.size());
ParallelChunks(n_pair, ThreadsForWork(anom.size(), nthreads, 4096), [&](int lo, int hi) {
for (int q = lo; q < hi; ++q)
pair_centric[q] = gops.is_reflection_centric(
gemmi::Op::Miller{pair_hkl[q].h, pair_hkl[q].k, pair_hkl[q].l}) ? 1 : 0;
});
for (size_t p = 0; p < anom.size(); ++p) {
AnomAcc a = anom[p];
a.h = pair_hkl[p].h; a.k = pair_hkl[p].k; a.l = pair_hkl[p].l;
const uint64_t fkey = pair_hkl[p].pack();
const bool centric = gops.is_reflection_centric(gemmi::Op::Miller{a.h, a.k, a.l});
const bool centric = pair_centric[p] != 0;
{
const auto hshell = shells.GetShell(a.d);
if (hshell && *hshell >= 0 && *hshell < n_shells) {
@@ -6326,7 +6433,10 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
restored_on_gpu = true;
}
#endif
if (!restored_on_gpu)
// The CPU partial scaling below writes every partial's corr when it is done, from a copy it starts
// off corr_ingested itself, and nothing between here and there reads corr - so there the restore
// is not made, which saves a pass over the eighty-byte records to set four bytes of each.
if (!restored_on_gpu && !scale_partials)
ParallelChunks(static_cast<int>(partials.size()), ThreadsForWork(partials.size(), nthreads),
[&](int lo, int hi) {
for (int i = lo; i < hi; ++i) partials[i].corr = corr_ingested[i];
@@ -6416,21 +6526,51 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
} else {
// The loop works on its own compact copy of the partials (ScalingRef and ScalingFit in the
// header) and hands corr back to them when it is done. scaling_index[q] is the partial at
// position q of the copy: the CSR's group order, then the partials in no group.
// position q of the copy: first the partials a group mean can take in - in a group and not
// below the partiality floor (ReferenceWeight's fixed tests) - group by group in the CSR's
// order, then all the others in index order. The means walk only the first part: a partial
// they would skip in any case adds nothing to a sum, so leaving it out of the walk changes no
// sum, and on a fine-sliced sweep half of the partials are rocking-curve tails below the floor.
const int n_partials = static_cast<int>(partials.size());
const size_t partial_threads = ThreadsForWork(partials.size(), nthreads);
const int n_grouped = partial_group_start.empty()
? 0 : partial_group_start.back() + partial_group_count.back();
std::vector<int32_t> scaling_index(partial_group_perm.get(), partial_group_perm.get() + n_grouped);
for (int i = 0; i < n_partials; ++i)
if (partials[i].group < 0) scaling_index.push_back(i);
std::vector<ScalingRef> scaling_refs(n_partials);
std::vector<float> scaling_corr(n_partials);
// Taken once, in index order, so that the walks over the CSR below gather one byte a partial
// rather than its record.
FillVector<uint8_t> in_ref(n_partials);
ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
for (int i = lo; i < hi; ++i)
in_ref[i] = partial_group[i] >= 0 && !(partials[i].partiality < min_partiality) ? 1 : 0;
});
auto in_reference = [&](int i) { return in_ref[i] != 0; };
std::vector<int32_t> ref_start(n_groups), ref_count(n_groups);
ParallelChunks(n_groups, ThreadsForWork(partials.size(), nthreads), [&](int glo, int ghi) {
for (int g = glo; g < ghi; ++g) {
int32_t c = 0;
for (int p = partial_group_start[g]; p < partial_group_start[g] + partial_group_count[g]; ++p)
c += in_reference(partial_group_perm[p]) ? 1 : 0;
ref_count[g] = c;
}
});
int n_ref = 0;
for (int g = 0; g < n_groups; ++g) { ref_start[g] = n_ref; n_ref += ref_count[g]; }
FillVector<int32_t> scaling_index(n_partials);
ParallelChunks(n_groups, ThreadsForWork(partials.size(), nthreads), [&](int glo, int ghi) {
for (int g = glo; g < ghi; ++g) {
int32_t q = ref_start[g];
for (int p = partial_group_start[g]; p < partial_group_start[g] + partial_group_count[g]; ++p)
if (in_reference(partial_group_perm[p])) scaling_index[q++] = partial_group_perm[p];
}
});
const std::vector<int32_t> rest = KeptIndices(n_partials, nthreads, [&](int i) { return !in_reference(i); });
ParallelChunks(static_cast<int>(rest.size()), ThreadsForWork(rest.size(), nthreads), [&](int lo, int hi) {
for (int j = lo; j < hi; ++j) scaling_index[n_ref + j] = rest[j];
});
FillVector<ScalingRef> scaling_refs(n_partials);
FillVector<float> scaling_corr(n_partials);
ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
for (int q = lo; q < hi; ++q) {
const Obs &o = partials[scaling_index[q]];
scaling_refs[q] = {o.frame, o.I, o.sigma, o.partiality, o.prescaling_corr};
scaling_corr[q] = o.corr;
scaling_corr[q] = corr_ingested[scaling_index[q]]; // the restore (see the top)
}
});
@@ -6439,68 +6579,74 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
// partials are rocking-curve tails below it, measured at I/sigma ~ 0.1, and a frame's scale
// fitted on them is fitted on noise that the reference does not answer to. The Rotation
// model needs zeta > 0.
auto fitted_obs = [&](const Obs &o) {
return o.group >= 0 && !o.on_ice && o.partiality >= min_partiality
auto fitted_obs = [&](int i) {
const Obs &o = partials[i];
return partial_group[i] >= 0 && !o.on_ice && o.partiality >= min_partiality
&& std::isfinite(o.zeta) && o.zeta > 0.0f;
};
const int n_fit_frames = static_cast<int>(frame_start.size());
std::vector<int32_t> fit_start(n_fit_frames + 1, 0);
ParallelFor(n_fit_frames, nthreads, [&](int f) {
fit_start[f + 1] = static_cast<int32_t>(std::count_if(
partials.begin() + frame_start[f], partials.begin() + frame_start[f] + frame_count[f],
fitted_obs));
int32_t n = 0;
for (int i = frame_start[f]; i < frame_start[f] + frame_count[f]; ++i) n += fitted_obs(i) ? 1 : 0;
fit_start[f + 1] = n;
});
for (int f = 0; f < n_fit_frames; ++f) fit_start[f + 1] += fit_start[f];
std::vector<ScalingFit> scaling_fit(fit_start.back());
FillVector<ScalingFit> scaling_fit(fit_start.back());
ParallelFor(n_fit_frames, nthreads, [&](int f) {
int32_t k = fit_start[f];
for (int i = frame_start[f]; i < frame_start[f] + frame_count[f]; ++i) {
if (!fitted_obs(i)) continue;
const Obs &o = partials[i];
if (!fitted_obs(o)) continue;
const double weight = SafeInv(o.sigma, 1.0);
scaling_fit[k++] = {o.group, o.partiality * SafeInv(o.prescaling_corr, 1.0),
static_cast<double>(o.I), weight * weight};
scaling_fit[k++] = {partial_group[i], o.I, o.sigma, o.partiality, o.prescaling_corr};
}
});
// RunScalingLoop rescales after every iteration before anything reads corr, so the iteration's
// UpdateCorr is folded into that pass over the partials: the G and the fitted frames as the
// iteration left them, then the rescale on top - the same two roundings, one pass instead of two.
std::vector<double> g_fitted;
std::vector<uint8_t> frame_fitted;
// UpdateCorr is folded into the rescale: the G and the fitted frames as the iteration left
// them, then the rescale on top - the same two roundings. Nor is that a pass of its own over
// the copy: the next iteration's group means are the next thing to read corr, so each partial
// they walk is brought up to date there, in the same walk, and whatever is still owed when the
// loop stops is applied after it. The partials the means do not walk are read by nothing until
// the loop is over, so their updates are kept - a few per-frame arrays an iteration - and
// applied to them once, in order, at the end. Each corr goes through the same updates in the
// same order.
CorrUpdate update;
bool update_owed = false; // the walked partials have not had `update` yet
std::vector<CorrUpdate> rest_updates;
partial_loop = RunScalingLoop("partials", g_partial, frame_count, smooth_window,
[&] {
ReduceScalingGroupMeans(n_groups, partial_group_start.data(), partial_group_count.data(),
scaling_refs, scaling_corr, partial_mean);
ReduceScalingGroupMeans(n_groups, ref_start.data(), ref_count.data(),
scaling_refs, scaling_corr, partial_mean, update_owed ? &update : nullptr);
update_owed = false;
FitScalingFrames(scaling_fit, fit_start, MIN_REFLECTIONS, partial_mean, g_partial);
g_fitted = g_partial;
frame_fitted = frame_scaled_scratch;
update.g = g_partial;
update.fitted = frame_scaled_scratch;
}, [&](const std::vector<uint8_t> &apply, const std::vector<double> &ratio) {
ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
for (int q = lo; q < hi; ++q) {
const ScalingRef &o = scaling_refs[q];
float &corr = scaling_corr[q];
if (frame_fitted[o.frame]) { // UpdateCorr
const double denom = static_cast<double>(o.partiality) * g_fitted[o.frame];
if (std::isfinite(o.prescaling_corr) && std::isfinite(denom) && denom > 0.0)
corr = static_cast<float>(o.prescaling_corr / denom);
else
corr = NAN;
}
if (apply[o.frame] && std::isfinite(corr))
corr = static_cast<float>(corr * ratio[o.frame]);
}
update.apply = apply;
update.ratio = ratio;
update_owed = true;
rest_updates.push_back(update);
});
if (update_owed)
ParallelChunks(n_ref, partial_threads, [&](int lo, int hi) {
for (int q = lo; q < hi; ++q) update.Apply(scaling_refs[q], scaling_corr[q]);
});
ParallelChunks(n_partials - n_ref, partial_threads, [&](int lo, int hi) {
for (int q = n_ref + lo; q < n_ref + hi; ++q)
for (const CorrUpdate &u : rest_updates) u.Apply(scaling_refs[q], scaling_corr[q]);
});
// Back into the partials' own order - the order the memo keeps - in the compact array, and
// from there onto the partials in one sequential pass.
FillVector<float> memo_corr(n_partials);
ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
for (int q = lo; q < hi; ++q) partials[scaling_index[q]].corr = scaling_corr[q];
for (int q = lo; q < hi; ++q) memo_corr[scaling_index[q]] = scaling_corr[q];
});
// The memo keeps corr in the partials' own order.
ParallelChunks(n_partials, partial_threads, [&](int lo, int hi) {
for (int i = lo; i < hi; ++i) scaling_corr[i] = partials[i].corr;
for (int i = lo; i < hi; ++i) partials[i].corr = memo_corr[i];
});
PartialLoopMemo m{key, partial_loop, frame_scaled_scratch, g_partial, std::move(scaling_corr)};
PartialLoopMemo m{key, partial_loop, frame_scaled_scratch, g_partial, std::move(memo_corr)};
if (partial_loop_memos.size() == MAX_PARTIAL_LOOP_MEMOS)
partial_loop_memos.erase(partial_loop_memos.begin());
partial_loop_memos.push_back(std::move(m));
@@ -6785,11 +6931,15 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
// full_index[q] is the full at position q, in the group CSR's order and then the fulls in no group.
const int nf = static_cast<int>(fulls.size());
const size_t full_threads = ThreadsForWork(fulls.size(), nthreads);
std::vector<int32_t> full_index(fulls_group_perm.begin(), fulls_group_perm.end());
for (int i = 0; i < nf; ++i)
if (fulls[i].group < 0) full_index.push_back(i);
std::vector<ScalingRef> full_refs(nf);
std::vector<float> full_corr(nf);
const int nf_grouped = static_cast<int>(fulls_group_perm.size());
const std::vector<int32_t> ungrouped = KeptIndices(nf, nthreads, [&](int i) { return fulls[i].group < 0; });
FillVector<int32_t> full_index(nf);
ParallelChunks(nf, full_threads, [&](int lo, int hi) {
for (int q = lo; q < hi; ++q)
full_index[q] = q < nf_grouped ? fulls_group_perm[q] : ungrouped[q - nf_grouped];
});
FillVector<ScalingRef> full_refs(nf);
FillVector<float> full_corr(nf);
ParallelChunks(nf, full_threads, [&](int lo, int hi) {
for (int q = lo; q < hi; ++q) {
const Obs &o = fulls[full_index[q]];
@@ -6812,14 +6962,13 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
fulls.begin() + fulls_frame_start[f] + fulls_frame_count[f], fitted_full));
});
for (int f = 0; f < n_frames; ++f) fit_start[f + 1] += fit_start[f];
std::vector<ScalingFit> full_fit(fit_start.back());
FillVector<ScalingFit> full_fit(fit_start.back());
ParallelFor(n_frames, nthreads, [&](int f) {
int32_t k = fit_start[f];
for (int i = fulls_frame_start[f]; i < fulls_frame_start[f] + fulls_frame_count[f]; ++i) {
const Obs &o = fulls[i];
if (!fitted_full(o)) continue;
const double weight = SafeInv(o.sigma, 1.0);
full_fit[k++] = {o.group, 1.0, static_cast<double>(o.I), weight * weight};
full_fit[k++] = {o.group, o.I, o.sigma, 1.0f, 1.0f};
}
});
image_analysis/scale_merge/RotationScaleMerge.h
+67 -19
View File
@@ -8,7 +8,11 @@
#include <cstdint>
#include <functional>
#include <limits>
#include <memory>
#include <new>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
#include "../../common/DiffractionExperiment.h"
@@ -26,6 +30,24 @@
#include "RotationScaleMergeGPU.h"
#endif
// A vector allocator that leaves each new element default-initialised - for a type without initialisers
// of its own, not written at all - instead of zeroing it. An array sized on one thread and then filled
// on all of them is otherwise first written page by page by that one thread: on a long sweep, gigabytes
// of page faults taken serially before the parallel fill starts. Everything else is std::allocator's.
template <class T>
struct DefaultInitAllocator : std::allocator<T> {
template <class U> struct rebind { using other = DefaultInitAllocator<U>; };
DefaultInitAllocator() noexcept = default;
template <class U> DefaultInitAllocator(const DefaultInitAllocator<U> &) noexcept {}
template <class U> void construct(U *p) noexcept(std::is_nothrow_default_constructible_v<U>) {
::new (static_cast<void *>(p)) U;
}
template <class U, class... Args> void construct(U *p, Args &&... args) {
::new (static_cast<void *>(p)) U(std::forward<Args>(args)...);
}
};
template <class T> using FillVector = std::vector<T, DefaultInitAllocator<T>>;
// Dedicated, allocate-once scale+combine+merge for rotation data (the -P rot3d path): recompute the
// per-frame partiality from the (smoothed) mosaicity, robustly fit a per-image scale G, 3D-combine each
// rocking event's partials into fulls, refit a per-frame scale on the fulls (XDS order), and merge with
@@ -230,17 +252,20 @@ private:
// One integrated observation - a per-frame partial during scaling/combine, or a combined full during
// scale-fulls/merge. Flat (not nested per image); a POD so the arrays translate straight to CUDA.
// No member initialisers, so that `partials` and `fulls` (FillVector) are first written by the
// parallel passes that fill them: whoever makes one sets every field.
struct Obs {
int32_t h, k, l;
float I, sigma, d, prescaling_corr, partiality, zeta, delta_phi, bkg, var_bkg;
// Fulls only, written by the combine: the full's variance as a function of intensity,
// var(I) = var_bkg + var_per_I * I. The merge rebuilds it at the reflection's mean.
float var_per_I = 0.0f;
// var(I) = var_bkg + var_per_I * I. The merge rebuilds it at the reflection's mean. 0 on a partial.
float var_per_I;
// Fulls only: the capture uncertainty, a relative sigma - capture_uncertainty_coeff times the
// part of the rocking curve the event did not catch - so the merge can add (capture * <I>)^2 to
// the variance it rebuilds at the reflection's mean, as the combine added it to the full's own.
float capture = 0.0f;
float px = NAN, py = NAN; // predicted detector position (for the absorption surface; CPU path only)
// 0 on a partial.
float capture;
float px, py; // predicted detector position (for the absorption surface; CPU path only)
float image_number; // fractional frame position (for 3D-combine contiguity)
int32_t frame; // index of the outcome whose per-frame scale G applies to this obs
uint8_t on_ice;
@@ -343,7 +368,7 @@ private:
// built only when the device pipeline is NOT resident (no GPU, or an observation dump was asked
// for); on the resident path (resident_ingest) the device holds the observations and the host
// keeps only the ingest arrays below, until the end of Ingest.
std::vector<Obs> partials; // all per-frame partials, grouped by frame
FillVector<Obs> partials; // all per-frame partials, grouped by frame
bool resident_ingest = false;
// What the ingest-time geometry smoothing (SmoothMosaicityAndPartiality) needs, one array per field
// in `partials` order: the two fields it rewrites, the two the rocking-event walks read in raw-hkl
@@ -420,7 +445,7 @@ private:
std::vector<float> rawrun_d; // representative resolution per raw hkl
std::vector<int32_t> rawrun_group; // dense ASU-group id per raw hkl (<0 = absent/out of range)
std::vector<Obs> fulls; // combined fulls (rebuilt each Run), sorted by frame
FillVector<Obs> fulls; // combined fulls (rebuilt each Run), sorted by frame
// The fulls by ASU group, each group's in index order: group g is fulls_group_perm[fulls_group_start[g]]
// up to fulls_group_start[g + 1]. A pass over every group's fulls can then give each thread whole
// groups and still add a group up in the order a serial walk over the fulls would. Rebuilt with the
@@ -624,6 +649,10 @@ private:
// group over it, on all threads and in the order the serial loop would.
std::vector<int32_t> partial_group_start, partial_group_count;
NotCopied<int32_t[]> partial_group_perm;
// Without a GPU: each partial's ASU group (-1 = in none), in observation order - where the CPU
// passes read a partial's group from. A partial's own Obs::group is not stamped: written into the
// eighty-byte records at random, once per grouping, it cost more than the grouping itself.
NotCopied<int32_t[]> partial_group;
// The CPU scaling loops' working copy of the partials, or of the fulls. A loop passes over them
// several times per iteration and reads a handful of fields each time, so it works on these rather
@@ -637,10 +666,13 @@ private:
int32_t frame;
float I, sigma, partiality, prescaling_corr;
};
// The fit's factors that do not depend on the reference - the coefficient of a mean of 1,
// partiality / prescaling_corr (1 for a full), and the weight 1/sigma^2 - are worked out from these
// floats where they are used: twenty bytes an observation read every iteration instead of thirty-two,
// and the same doubles.
struct ScalingFit {
int32_t group;
double coeff_per_mean; // partiality / prescaling_corr (1 for a full): the fit's coeff for a mean of 1
double Iobs, w2;
float I, sigma, partiality, prescaling_corr;
};
// The observation's inverse-variance weight in the scaling reference, and its scaled intensity;
@@ -656,17 +688,33 @@ private:
I_corr = Ic;
return 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
}
static double ReferenceWeight(const Obs &o, double min_partiality, float &I_corr) {
return ReferenceWeight(o.group, o.corr, o.partiality, o.I, o.sigma, min_partiality, I_corr);
}
// Inverse-variance per-group mean of I*corr over the partials (the scaling reference; no cell mask -
// only the final merge applies one), using the CSR above, each group summed in index order.
void ReducePartialGroupMeans(int n_groups, std::vector<double> &out_mean) const;
// The partials' loop's corr update of one iteration - UpdateCorr from the G it fitted, then the
// rescale RunScalingLoop asked for - per frame: `fitted` and `g` as the iteration left them, `apply`
// and `ratio` the rescale's.
struct CorrUpdate {
std::vector<uint8_t> fitted, apply;
std::vector<double> g, ratio;
void Apply(const ScalingRef &o, float &corr) const {
if (fitted[o.frame]) { // UpdateCorr
const double denom = static_cast<double>(o.partiality) * g[o.frame];
if (std::isfinite(o.prescaling_corr) && std::isfinite(denom) && denom > 0.0)
corr = static_cast<float>(o.prescaling_corr / denom);
else
corr = NAN;
}
if (apply[o.frame] && std::isfinite(corr))
corr = static_cast<float>(corr * ratio[o.frame]);
}
};
// The same again over a scaling loop's copy, refs and corr both in group order: group g's at
// [group_start[g], group_start[g] + group_count[g]).
// [group_start[g], group_start[g] + group_count[g]). With `update`, each corr is first brought up to
// date by it, in the same walk (see the partials' loop in RunArm).
void ReduceScalingGroupMeans(int n_groups, const int32_t *group_start, const int32_t *group_count,
const std::vector<ScalingRef> &refs, const std::vector<float> &corr,
std::vector<double> &out_mean) const;
const FillVector<ScalingRef> &refs, FillVector<float> &corr,
std::vector<double> &out_mean, const CorrUpdate *update = nullptr) const;
// Per-frame G fit against the reference group_mean_in[group], over a scaling loop's ScalingFit, frame
// f's at [fit_start[f], fit_start[f + 1]). A frame with fewer than min_obs of them keeps its G: a
@@ -674,7 +722,7 @@ private:
// event filed under one frame, so on a sparse sweep a frame holds a handful - its fit goes in at
// the information it has (min_obs 1), and the smoothing in RunFullsScalingLoop decides how far its
// neighbours carry it.
void FitScalingFrames(const std::vector<ScalingFit> &fit, const std::vector<int32_t> &fit_start,
void FitScalingFrames(const FillVector<ScalingFit> &fit, const std::vector<int32_t> &fit_start,
size_t min_obs, const std::vector<double> &group_mean_in, std::vector<double> &g);
// The alternating per-frame scale fit run to convergence: `iterate` fits every frame once against
@@ -710,7 +758,7 @@ private:
ScalingLoopOutcome outcome;
std::vector<uint8_t> fitted;
std::vector<double> g;
std::vector<float> corr;
FillVector<float> corr;
};
// A pass alternates its groupings (P1 full range, P1 at the search limit, the adopted group), so a
// few are kept; the oldest go first.
@@ -863,9 +911,9 @@ private:
const char *name, bool resolution_gauge,
const SphericalHarmonicBasis *sh_basis = nullptr);
// Sort `fulls` by peak frame and (re)build fulls_frame_start/count (the per-frame CSR the scale-fulls
// step slices). Shared by the CPU Combine tail and the GPU combine path.
void SortFullsByFrame();
// Fill `fulls` with the parts' fulls sorted by peak frame and (re)build fulls_frame_start/count (the
// per-frame CSR the scale-fulls step slices). The CPU Combine's tail.
void SortFullsByFrame(const std::vector<std::vector<Obs>> &part);
// Per-frame CC vs the partial merge reference (CPU; the GPU equivalent is gpu_->ComputePartialCC).
void ComputePerFrameCC(const std::vector<double> &partial_group_mean,