feat(data): SPSC ring buffer, incremental min/max LOD pyramid, series types

src/data/minmax_pyramid.cpp

@@ -0,0 +1,240 @@
+#include <plot/data/minmax_pyramid.hpp>
+
+#include <algorithm>
+#include <cmath>
+
+namespace plot::data {
+namespace {
+
+/// Samples covered by one level-l bucket.
+constexpr u64 level_span(u32 l) { return u64{kBucket} << l; }
+
+/// Level-l buckets needed to hold n samples (ceiling division).
+constexpr u64 buckets_for(u64 n, u32 l) {
+  const u64 span = level_span(l);
+  return n / span + (n % span == 0 ? u64{0} : u64{1});
+}
+
+/// First level-0 bucket boundary at or after v.
+constexpr u64 ceil_to_bucket(u64 v) {
+  return (v / kBucket + (v % kBucket == 0 ? u64{0} : u64{1})) * kBucket;
+}
+
+/// Canonical, coarsest-first decomposition of the level-0 bucket span [b0, b1)
+/// into pyramid entries, calling fn(level, index) for each.
+///
+/// A level-l entry is emitted only when the whole level-l bucket exists (all
+/// kBucket << l samples are present, so its aggregate is final). Callers
+/// guarantee b1 <= ceil(size() / kBucket), so the only possible remainder is the
+/// single trailing partial level-0 bucket, which is emitted at level 0. The
+/// decomposition has O(log) entries, so every query below is O(log n).
+template <typename Fn>
+void cover_buckets(const MinMaxPyramid& p, u64 b0, u64 b1, Fn&& fn) {
+  const u64 complete = p.size() / kBucket;
+  u64 b = b0;
+  for (u32 li = p.level_count(); li-- > 0 && b < b1;) {
+    const u64 lb = u64{1} << li;               // level-0 buckets per level-li bucket
+    const u64 avail = p.level(li).size();
+    while (b < b1 && (b % lb) == 0 && b + lb <= b1 && b + lb <= complete &&
+           (b / lb) < avail) {
+      fn(li, b / lb);
+      b += lb;
+    }
+  }
+  if (b < b1) fn(0u, b);                       // trailing partial level-0 bucket
+}
+
+/// Folds the decomposition of [b0, b1) into a min/max Range. Returns a
+/// zero-span Range when there is nothing to fold.
+Range fold_range(const MinMaxPyramid& p, u64 b0, u64 b1) {
+  Range r{};
+  bool first = true;
+  cover_buckets(p, b0, b1, [&](u32 l, u64 i) {
+    const MinMax m = p.level(l)[static_cast<usize>(i)];
+    if (first) {
+      r.min = f64(m.min);
+      r.max = f64(m.max);
+      first = false;
+    } else {
+      r.min = std::min(r.min, f64(m.min));
+      r.max = std::max(r.max, f64(m.max));
+    }
+  });
+  return first ? Range{0.0, 0.0} : r;
+}
+
+}  // namespace
+
+void MinMaxPyramid::add_level() {
+  levels_.emplace_back();
+  sums_.emplace_back();
+  sumsqs_.emplace_back();
+}
+
+void MinMaxPyramid::reserve(u64 n) {
+  if (n > reserved_n_) reserved_n_ = n;
+  const u64 want = std::max(std::max(n, size_), u64{2});
+
+  // Materialise every level n samples can need -- buckets *and* prefix sums --
+  // and size each one up front, so a run of appends allocates nothing at all.
+  u32 max_levels = 1;
+  for (u64 span = level_span(0); span <= want && max_levels < 63; span <<= 1) ++max_levels;
+  while (levels_.size() < max_levels) add_level();
+  for (u32 l = 0; l < levels_.size(); ++l) {
+    const usize want_l = static_cast<usize>(buckets_for(want, l));
+    levels_[l].reserve(want_l);
+    sums_[l].reserve(want_l);
+    sumsqs_[l].reserve(want_l);
+  }
+}
+
+void MinMaxPyramid::append(f32 v) {
+  const u64 prev = size_;
+  const f64 vd = f64(v);
+  // sum_[l].back() is the sum over every sample appended so far, and the newest
+  // sample always lands in the last bucket of every level, so the last prefix
+  // entry of every level is exactly the running total.
+  const f64 total_before = sums_[0].empty() ? 0.0 : sums_[0].back();
+  const f64 total = total_before + vd;
+  const f64 sq_before = sumsqs_[0].empty() ? 0.0 : sumsqs_[0].back();
+  const f64 total_sq = sq_before + vd * vd;
+
+  const u64 idx = prev / kBucket;
+  if (idx == levels_[0].size()) {
+    levels_[0].push_back(MinMax{v, v});
+    sums_[0].push_back(total);
+    sumsqs_[0].push_back(total_sq);
+  } else {
+    MinMax& cur = levels_[0][static_cast<usize>(idx)];
+    cur.min = std::min(cur.min, v);
+    cur.max = std::max(cur.max, v);
+    sums_[0][static_cast<usize>(idx)] = total;
+    sumsqs_[0][static_cast<usize>(idx)] = total_sq;
+  }
+
+  for (u32 l = 1; l < 63; ++l) {
+    // A level only carries data once the level below has more than one bucket: a
+    // one-child parent would just duplicate it.
+    if (levels_[l - 1].size() < 2) break;
+    if (l == live_levels_) {
+      ++live_levels_;
+      if (l == levels_.size()) {   // no reserve() sized this level up front
+        add_level();
+        if (reserved_n_ != 0) {
+          const usize want = static_cast<usize>(buckets_for(reserved_n_, l));
+          levels_[l].reserve(want);
+          sums_[l].reserve(want);
+          sumsqs_[l].reserve(want);
+        }
+      }
+    }
+    const u64 ci = prev / level_span(l - 1);   // current bucket at level l-1
+    const u64 parent = ci / 2;                 // its parent at level l
+    MinMax m;
+    if ((ci & 1u) != 0) {
+      const MinMax& a = levels_[l - 1][static_cast<usize>(ci - 1)];
+      const MinMax& b = levels_[l - 1][static_cast<usize>(ci)];
+      m.min = std::min(a.min, b.min);
+      m.max = std::max(a.max, b.max);
+    } else {
+      m = levels_[l - 1][static_cast<usize>(ci)];
+    }
+    if (parent == levels_[l].size()) {
+      levels_[l].push_back(m);
+      sums_[l].push_back(total);
+      sumsqs_[l].push_back(total_sq);
+    } else {
+      levels_[l][static_cast<usize>(parent)] = m;
+      sums_[l][static_cast<usize>(parent)] = total;
+      sumsqs_[l][static_cast<usize>(parent)] = total_sq;
+    }
+  }
+  size_ = prev + 1;
+}
+
+void MinMaxPyramid::clear() {
+  size_ = 0;
+  live_levels_ = 1;
+  for (auto& lv : levels_) lv.clear();
+  for (auto& s : sums_) s.clear();
+  for (auto& s : sumsqs_) s.clear();
+}
+
+std::span<const MinMax> MinMaxPyramid::level(u32 l) const {
+  if (l >= level_count()) return {};
+  return levels_[l];
+}
+
+Range MinMaxPyramid::reduce_buckets(u64 begin, u64 end) const {
+  const u64 e = std::min(end, size_);
+  const u64 b = std::min(begin, e);
+  const u64 b0 = ceil_to_bucket(b) / kBucket;   // first whole bucket inside
+  const u64 b1 = e / kBucket;                   // one past the last whole bucket
+  if (b1 <= b0) return Range{0.0, 0.0};
+  return fold_range(*this, b0, b1);
+}
+
+f64 MinMaxPyramid::sum_buckets(u64 begin, u64 end) const {
+  const u64 e = std::min(end, size_);
+  const u64 b = std::min(begin, e);
+  const u64 b0 = ceil_to_bucket(b) / kBucket;
+  const u64 b1 = e / kBucket;
+  if (b1 <= b0) return 0.0;
+  f64 total = 0.0;
+  cover_buckets(*this, b0, b1, [&](u32 l, u64 i) {
+    const f64 prev = i == 0 ? 0.0 : sums_[l][static_cast<usize>(i - 1)];
+    total += sums_[l][static_cast<usize>(i)] - prev;
+  });
+  return total;
+}
+
+f64 MinMaxPyramid::sumsq_buckets(u64 begin, u64 end) const {
+  const u64 e = std::min(end, size_);
+  const u64 b = std::min(begin, e);
+  const u64 b0 = ceil_to_bucket(b) / kBucket;
+  const u64 b1 = e / kBucket;
+  if (b1 <= b0) return 0.0;
+  f64 total = 0.0;
+  cover_buckets(*this, b0, b1, [&](u32 l, u64 i) {
+    const f64 prev = i == 0 ? 0.0 : sumsqs_[l][static_cast<usize>(i - 1)];
+    total += sumsqs_[l][static_cast<usize>(i)] - prev;
+  });
+  return total;
+}
+
+Range MinMaxPyramid::reduce(u64 begin, u64 end) const {
+  const u64 e = std::min(end, size_);
+  const u64 b = std::min(begin, e);
+  if (e <= b) return Range{0.0, 0.0};
+  const u64 b0 = b / kBucket;                   // first bucket holding a sample
+  const u64 b1 = (e - 1) / kBucket + 1;         // one past the last one
+  return fold_range(*this, b0, b1);
+}
+
+u32 MinMaxPyramid::level_for(u64 n, f64 samples_per_pixel) const {
+  // The ratio test also rejects NaN, and the clamp keeps a huge (or infinite)
+  // samples-per-pixel from making the floor-to-u32 conversion ill-formed.
+  const f64 ratio = samples_per_pixel / f64(kBucket);
+  if (n == 0 || !(ratio > 1.0)) return 0;
+  const f64 lg = std::floor(std::log2(ratio));
+  u32 l = lg >= f64(62) ? 62u : static_cast<u32>(lg);
+  while (l > 0) {
+    const bool exists = l < level_count() && level(l).size() >= 2;
+    if (exists && buckets_for(n, l) >= 2) break;
+    --l;
+  }
+  return l;
+}
+
+u64 MinMaxPyramid::snap(f64 x, f64 x0, f64 dx) const {
+  if (dx <= 0.0 || size_ == 0) return 0;
+  const f64 i = (x - x0) / dx;
+  // Bound the round before it: llround() is undefined for values a long long
+  // cannot hold, and this also maps NaN to index 0.
+  if (!(i > 0.0)) return 0;
+  if (i >= f64(size_)) return size_ - 1;
+  const u64 idx = static_cast<u64>(std::llround(i));
+  return idx >= size_ ? size_ - 1 : idx;
+}
+
+}  // namespace plot::data