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

WilsonOutliers: judged list, per-reflection verdicts and both sorts on all threads

Browse Source 
Exact (tier E): the same observations in the same order, the same verdicts.

- The judged observations are compacted in fixed blocks (count, prefix, fill) instead of a
  serial push_back.
- Both sorts (by d, by reflection) sort (key, index) records, so a comparison reads the two
  elements it compares instead of two observations elsewhere in memory; same total order.
- Each reflection is decided on its own observations, so the reflections go side by side in
  claimed blocks; the rejection count is an integer sum.

Measured (prototype, GPU): Wilson test 8tyy 3.4 -> 3.0 s, 8a1a 2.8 -> 2.6 s per run (the
sort-on-records part is not in the measured binary yet); md5 identical.

Clean branch rebuilt from scratch (GPU and CPU) and re-verified: p.mtz md5 identical to production on myob/cytc/thau/8a1a/8tyy/8qaw/9gdj (GPU), cytc -N 4, myob/cytc and myob -N 4 (CPU).

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 1f5c63f13c
commit 70859bcb60
1 changed files with 91 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/WilsonOutliers.cpp
+91 -38
View File
@@ -14,13 +14,29 @@ WilsonOutlierResult WilsonOutliers(const std::vector<WilsonObservation> &obs, do
out.rejected.assign(obs.size(), 0);
out.e2.assign(obs.size(), NAN);
std::vector<size_t> idx;
for (size_t i = 0; i < obs.size(); ++i) {
// The observations that can be judged, in index order: fixed blocks count theirs, a prefix places
// them, and every block writes its own stretch.
const auto judged = [&](size_t i) {
const auto &o = obs[i];
if (o.unit >= 0 && o.d > 0.0f && o.d < WILSON_OUTLIER_D_MAX && std::isfinite(o.I)
&& o.sigma > 0.0f && std::isfinite(o.sigma))
idx.push_back(i);
}
return o.unit >= 0 && o.d > 0.0f && o.d < WILSON_OUTLIER_D_MAX && std::isfinite(o.I)
&& o.sigma > 0.0f && std::isfinite(o.sigma);
};
const int n_obs = static_cast<int>(obs.size());
const int nb = ReductionBlocks(n_obs, 16384);
const auto block_begin = [&](int b) { return static_cast<size_t>(static_cast<int64_t>(n_obs) * b / nb); };
std::vector<size_t> block_at(nb + 1, 0);
ParallelFor(nb, nthreads, [&](int b) {
size_t c = 0;
for (size_t i = block_begin(b); i < block_begin(b + 1); ++i) c += judged(i) ? 1 : 0;
block_at[b + 1] = c;
});
for (int b = 0; b < nb; ++b) block_at[b + 1] += block_at[b];
std::vector<size_t> idx(nb > 0 ? block_at[nb] : 0);
ParallelFor(nb, nthreads, [&](int b) {
size_t k = block_at[b];
for (size_t i = block_begin(b); i < block_begin(b + 1); ++i)
if (judged(i)) idx[k++] = i;
});
if (idx.empty())
return out;
@@ -37,9 +53,21 @@ WilsonOutlierResult WilsonOutliers(const std::vector<WilsonObservation> &obs, do
// put the shell mean at 2% (an exponential's sd equals its mean); narrow shells keep the fall-off
// across one from reading as spread.
constexpr size_t OBS_PER_SHELL = 2000;
ParallelSort(idx.begin(), idx.end(), nthreads, [&](size_t a, size_t b) {
return obs[a].d != obs[b].d ? obs[a].d > obs[b].d : a < b;
});
// Sorted with the key beside the index, so a comparison reads the two elements it compares rather
// than two observations somewhere else in memory; the same order.
{
struct ByD { float d; size_t i; };
std::vector<ByD> by_d(idx.size());
ParallelChunks(static_cast<int>(idx.size()), nthreads, [&](int lo, int hi) {
for (int j = lo; j < hi; ++j) by_d[j] = {obs[idx[j]].d, idx[j]};
});
ParallelSort(by_d.begin(), by_d.end(), nthreads, [](const ByD &a, const ByD &b) {
return a.d != b.d ? a.d > b.d : a.i < b.i;
});
ParallelChunks(static_cast<int>(idx.size()), nthreads, [&](int lo, int hi) {
for (int j = lo; j < hi; ++j) idx[j] = by_d[j].i;
});
}
const size_t n_shells = std::max<size_t>(1, idx.size() / OBS_PER_SHELL);
// The observations in that order, so the walks below read them one after another instead of each
// through idx, and where each shell starts: shell s is the j with j * n_shells / N == s, which begins
@@ -144,41 +172,66 @@ WilsonOutlierResult WilsonOutliers(const std::vector<WilsonObservation> &obs, do
if (shell_valid[sh])
tested.insert(tested.end(), idx.begin() + shell_start[sh], idx.begin() + shell_start[sh + 1]);
out.n_tested = tested.size();
ParallelSort(tested.begin(), tested.end(), nthreads, [&](size_t a, size_t b) {
return obs[a].unit != obs[b].unit ? obs[a].unit < obs[b].unit : a < b;
});
{
struct ByUnit { int32_t unit; size_t i; };
std::vector<ByUnit> by_unit(tested.size());
ParallelChunks(static_cast<int>(tested.size()), nthreads, [&](int lo, int hi) {
for (int j = lo; j < hi; ++j) by_unit[j] = {obs[tested[j]].unit, tested[j]};
});
ParallelSort(by_unit.begin(), by_unit.end(), nthreads, [](const ByUnit &a, const ByUnit &b) {
return a.unit != b.unit ? a.unit < b.unit : a.i < b.i;
});
ParallelChunks(static_cast<int>(tested.size()), nthreads, [&](int lo, int hi) {
for (int j = lo; j < hi; ++j) tested[j] = by_unit[j].i;
});
}
// An improbable observation measured once is rejected. One with company is rejected when most of
// its reflection's other observations are probable, unclipped witnesses and it disagrees with their
// mean beyond the errors - for a pair, the other member. Where most are large, the reflection is
// (Aimless's "keep if most observations are large").
for (size_t lo = 0; lo < tested.size();) {
size_t hi = lo;
while (hi < tested.size() && obs[tested[hi]].unit == obs[tested[lo]].unit) ++hi;
for (size_t q = lo; q < hi; ++q) {
const size_t i = tested[q];
if (!improbable[i]) continue;
size_t n_probable = 0;
double sw = 0.0, swI = 0.0;
for (size_t r = lo; r < hi; ++r) {
const auto &o = obs[tested[r]];
if (r == q || improbable[tested[r]] || o.clipped) continue;
const double w = 1.0 / (static_cast<double>(o.sigma) * o.sigma);
sw += w; swI += w * o.I;
++n_probable;
}
const size_t n_others = hi - lo - 1;
bool reject = n_others == 0;
if (2 * n_probable > n_others) {
const auto &o = obs[i];
reject = std::fabs(o.I - swI / sw) > z * std::sqrt(static_cast<double>(o.sigma) * o.sigma + 1.0 / sw);
}
if (reject) {
out.rejected[i] = 1;
++out.n_rejected;
// Each reflection is decided on its own observations alone, so the reflections go side by side:
// where each one's stretch of `tested` starts, then the stretches in blocks.
const int n_tested = static_cast<int>(tested.size());
std::vector<uint8_t> starts_unit(n_tested);
ParallelChunks(n_tested, ThreadsForWork(tested.size(), nthreads), [&](int lo, int hi) {
for (int q = lo; q < hi; ++q)
starts_unit[q] = q == 0 || obs[tested[q]].unit != obs[tested[q - 1]].unit;
});
std::vector<int32_t> unit_start;
for (int q = 0; q < n_tested; ++q)
if (starts_unit[q]) unit_start.push_back(q);
const int n_units = static_cast<int>(unit_start.size());
unit_start.push_back(n_tested);
std::vector<size_t> rejected_in_block(ReductionBlocks(n_units, 1024), 0);
ParallelBlocks(n_units, nthreads, [&](int blk, int u0, int u1) {
for (int u = u0; u < u1; ++u) {
const size_t lo = unit_start[u], hi = unit_start[u + 1];
for (size_t q = lo; q < hi; ++q) {
const size_t i = tested[q];
if (!improbable[i]) continue;
size_t n_probable = 0;
double sw = 0.0, swI = 0.0;
for (size_t r = lo; r < hi; ++r) {
const auto &o = obs[tested[r]];
if (r == q || improbable[tested[r]] || o.clipped) continue;
const double w = 1.0 / (static_cast<double>(o.sigma) * o.sigma);
sw += w; swI += w * o.I;
++n_probable;
}
const size_t n_others = hi - lo - 1;
bool reject = n_others == 0;
if (2 * n_probable > n_others) {
const auto &o = obs[i];
reject = std::fabs(o.I - swI / sw) > z * std::sqrt(static_cast<double>(o.sigma) * o.sigma + 1.0 / sw);
}
if (reject) {
out.rejected[i] = 1;
++rejected_in_block[blk];
}
}
}
lo = hi;
}
}, 1024);
for (const size_t r : rejected_in_block) out.n_rejected += r;
return out;
}