diff --git a/image_analysis/scale_merge/WilsonOutliers.cpp b/image_analysis/scale_merge/WilsonOutliers.cpp index cfaaf8117..ef2c9e83c 100644 --- a/image_analysis/scale_merge/WilsonOutliers.cpp +++ b/image_analysis/scale_merge/WilsonOutliers.cpp @@ -14,13 +14,29 @@ WilsonOutlierResult WilsonOutliers(const std::vector &obs, do out.rejected.assign(obs.size(), 0); out.e2.assign(obs.size(), NAN); - std::vector 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(obs.size()); + const int nb = ReductionBlocks(n_obs, 16384); + const auto block_begin = [&](int b) { return static_cast(static_cast(n_obs) * b / nb); }; + std::vector 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 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 &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 by_d(idx.size()); + ParallelChunks(static_cast(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(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(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 &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 by_unit(tested.size()); + ParallelChunks(static_cast(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(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(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(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(tested.size()); + std::vector 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 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(unit_start.size()); + unit_start.push_back(n_tested); + std::vector 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(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(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; }