diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 94a9325af..aaaa11d02 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -7133,11 +7133,39 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b [](const auto &a, const auto &b) { return a.first > b.first; }); // low -> high res const int bins = std::clamp(static_cast(rs.size() / 100), 1, 40); const size_t per = (rs.size() + bins - 1) / static_cast(bins); - for (size_t b = 0; b * per < rs.size(); ++b) { - const size_t lo = b * per, hi = std::min(rs.size(), lo + per); + // The shell means, then their best NON-INCREASING fit (pool adjacent violators, weighted + // by shell size). The cut is where the FIT drops under 1, not where the first shell + // does: on a merge whose is flat near 1 - its ISa has collapsed - the first + // single shell to dip under 1 is a coin toss that moved the cut from 0.85 to 1.19 A + // between two compiler flags on one crystal, and the coarser cut left its glide zones + // too few absences to be judged. The fit moves by as much as its input does. + std::vector mean, weight; + for (size_t lo = 0; lo < rs.size(); lo += per) { + const size_t hi = std::min(rs.size(), lo + per); double sum = 0.0; for (size_t j = lo; j < hi; ++j) sum += rs[j].second; - if (sum / static_cast(hi - lo) < 1.0) { + mean.push_back(sum / static_cast(hi - lo)); + weight.push_back(static_cast(hi - lo)); + } + std::vector block_mean, block_weight; + std::vector block_size; + for (size_t b = 0; b < mean.size(); ++b) { + block_mean.push_back(mean[b]); block_weight.push_back(weight[b]); block_size.push_back(1); + while (block_mean.size() > 1 && block_mean[block_mean.size() - 2] < block_mean.back()) { + const size_t k = block_mean.size() - 1; + const double w = block_weight[k - 1] + block_weight[k]; + block_mean[k - 1] = (block_mean[k - 1] * block_weight[k - 1] + block_mean[k] * block_weight[k]) / w; + block_weight[k - 1] = w; + block_size[k - 1] += block_size[k]; + block_mean.pop_back(); block_weight.pop_back(); block_size.pop_back(); + } + } + std::vector fit; + for (size_t k = 0; k < block_mean.size(); ++k) + fit.insert(fit.end(), block_size[k], block_mean[k]); + for (size_t b = 0; b < fit.size(); ++b) { + const size_t lo = b * per; + if (fit[b] < 1.0) { // Cut the noise-dominated high-res shells. When even the lowest-res shell fails, // keep that shell alone rather than abandoning the cut. The bound is absolute // while the merged I/sigma it tests saturates at the merge's own ISa (times the