// SPDX-License-Identifier: MIT
// SPDX-FileCopyrightText: 2026 Martino Ferrari

#include "harness.hpp"
#include <plot/data/array_series.hpp>
#include <plot/data/detail/merge_count.hpp>
#include <plot/data/stream_series.hpp>
#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <random>
#include <thread>
#include <vector>

using plot::data::kBucket;

PLOT_TEST(array_series_bounds_and_gather) {
  std::vector<f32> y{1, 5, -2, 3};
  plot::data::ArraySeries s(std::move(y), /*x0=*/10.0, /*dx=*/0.5);
  PLOT_CHECK(s.size() == 4);
  PLOT_CHECK_NEAR(s.x_bounds().min, 10.0, 1e-9);
  PLOT_CHECK_NEAR(s.x_bounds().max, 11.5, 1e-9);
  PLOT_CHECK_NEAR(s.y_bounds().min, -2.0, 1e-6);
  std::vector<f32> out;
  s.gather(0, 0, 4, out);
  PLOT_CHECK(out.size() == 8);             // interleaved x,y
  // x is relative to x0(): sample 0 emits 0, not 10 (the relative-x contract).
  PLOT_CHECK_NEAR(out[0], 0.0, 1e-6);
  PLOT_CHECK_NEAR(out[1], 1.0, 1e-6);
  PLOT_CHECK_NEAR(out[6], 1.5, 1e-6);      // sample 3: x_at(3) - x0() = 3 * 0.5
}

PLOT_TEST(stream_series_acquire_publishes_snapshot) {
  plot::data::StreamSeries s(1 << 12, /*x0=*/0.0, /*dx=*/1e-6);
  for (int i = 0; i < 100; ++i) s.push(f32(i));
  PLOT_CHECK(s.produced() == 100);
  s.acquire();
  PLOT_CHECK(s.size() == 100);
  // Pyramid must track the published samples.
  PLOT_CHECK_NEAR(s.reduce(0, 100).max, 99.0, 1e-4);
  for (int i = 100; i < 250; ++i) s.push(f32(i));
  s.acquire();
  PLOT_CHECK(s.size() == 250);
}

PLOT_TEST(stream_series_drops_oldest_when_full) {
  plot::data::StreamSeries s(256, 0.0, 1.0);
  for (int i = 0; i < 1000; ++i) s.push(f32(i));
  s.acquire();
  PLOT_CHECK(s.size() <= 256);
  PLOT_CHECK_NEAR(s.y_bounds().max, 999.0, 1e-4);
}

PLOT_TEST(array_series_exact_range_queries_match_bruteforce) {
  // The strong test for the exactness split: arbitrary, deliberately
  // unaligned ranges must match a brute-force scan of the raw samples.
  std::mt19937 rng(7);
  std::uniform_real_distribution<float> d(-5.0f, 5.0f);
  std::vector<f32> y(4000);
  for (auto& v : y) v = d(rng);
  plot::data::ArraySeries s(y, 0.0, 1.0);
  for (u64 a = 0; a < 4000; a += 173) {
    for (u64 b = a + 1; b <= 4000; b += 331) {
      float lo = y[a], hi = y[a];
      double sum = 0.0, sq = 0.0;
      for (u64 i = a; i < b; ++i) {
        lo = std::min(lo, y[i]); hi = std::max(hi, y[i]);
        sum += y[i]; sq += double(y[i]) * double(y[i]);
      }
      const plot::Range r = s.reduce(a, b);
      PLOT_CHECK_NEAR(r.min, lo, 1e-4);
      PLOT_CHECK_NEAR(r.max, hi, 1e-4);
      PLOT_CHECK_NEAR(s.sum(a, b), sum, 1e-2);
      PLOT_CHECK_NEAR(s.sumsq(a, b), sq, 1e-1);
    }
  }
}

PLOT_TEST(stream_series_exact_queries_survive_ring_wrap) {
  plot::data::StreamSeries s(1u << 10, 0.0, 1.0);
  constexpr u64 kPushes = 3000;
  constexpr u64 kFirstRetained = kPushes - (1u << 10);  // 1976
  for (u64 i = 0; i < kPushes; ++i) s.push(f32(i % 97));   // wraps the ring twice
  s.acquire();
  const u64 n = s.size();
  PLOT_CHECK(n == 1u << 10);
  PLOT_CHECK_NEAR(s.x0(), 1976.0, 1e-9);
  PLOT_CHECK(s.x_bounds().min == s.x0());
  PLOT_CHECK(s.pyramid().size() == n);

  const plot::Range r = s.reduce(0, n);                 // retained window, wrapped
  PLOT_CHECK_NEAR(r.max, 96.0, 1e-4);
  PLOT_CHECK_NEAR(r.min, 0.0, 1e-4);

  double sum = 0.0, sumsq = 0.0;
  for (u64 i = kFirstRetained; i < kPushes; ++i) {
    const double v = static_cast<double>(i % 97);
    sum += v;
    sumsq += v * v;
  }
  PLOT_CHECK_NEAR(s.sum(0, n), sum, 1e-2);
  PLOT_CHECK_NEAR(s.sumsq(0, n), sumsq, 1e-1);

  std::vector<f32> out;
  s.gather(0, 0, n, out);                               // level 0: raw samples
  PLOT_CHECK(out.size() == 2 * n);
  // x is relative to x0(): the oldest retained sample emits 0, the newest
  // (n - 1) * dx, never the absolute data coordinates (relative-x contract).
  PLOT_CHECK_NEAR(out[0], 0.0, 1e-6);
  PLOT_CHECK_NEAR(out[1], f32(kFirstRetained % 97), 1e-6);
  PLOT_CHECK_NEAR(out[2 * n - 2], f64(n - 1) * 1.0, 1e-6);
}

PLOT_TEST(stream_series_push_and_read_are_race_free) {
  // The ring overwrites, so the producer reuses the slot of the snapshot's oldest
  // retained sample while the consumer may still be reading it. Relaxed atomic
  // slots make that defined; nothing else in the suite has two threads, so this is
  // the case ThreadSanitizer actually flags when the slots are plain.
  //
  // The two flags make the overlap structural, not scheduling luck: the reader
  // waits until the producer is already pushing, acknowledges, and only then does
  // the producer push the rest -- so the reader's window spans pushes that are in
  // flight by construction. The reader then loops until the producer reports done,
  // which also makes "at least one read" true by construction instead of by
  // winning a race with a fast producer.
  plot::data::StreamSeries s(1u << 8, 0.0, 1.0);   // one read spans reused slots
  constexpr u64 kStride = 512;                     // pushes before the reader starts
  constexpr u64 kPushes = 4'096;                   // total: a few thousand only
  std::atomic<bool> started{false};
  std::atomic<bool> reading{false};
  std::atomic<bool> done{false};

  std::thread producer([&] {
    for (u64 i = 0; i < kStride; ++i) s.push(f32(i % 971));
    started.store(true, std::memory_order_release);
    while (!reading.load(std::memory_order_acquire)) {
    }
    for (u64 i = kStride; i < kPushes; ++i) s.push(f32(i % 971));
    done.store(true, std::memory_order_release);
  });

  while (!started.load(std::memory_order_acquire)) {
  }
  reading.store(true, std::memory_order_release);

  u64 reads = 0;
  bool sizes_ok = true;
  bool finished = false;
  do {
    finished = done.load(std::memory_order_acquire);
    s.acquire();
    const u64 n = s.size();
    if (n > 0) {
      (void)s.reduce(0, n);
      (void)s.sum(0, n);
      (void)s.sumsq(0, n);
      std::vector<f32> out;
      s.gather(0, 0, n, out);                           // reads the oldest edge
      if (out.size() != 2 * n) sizes_ok = false;
    }
    ++reads;
  } while (!finished);

  producer.join();
  // Checked only after the join: a throw while the producer is joinable would
  // std::terminate instead of reporting the failure.
  PLOT_CHECK(s.produced() == kPushes);
  PLOT_CHECK(reads >= 1);
  PLOT_CHECK(sizes_ok);
}

namespace {

/// Buckets the pyramid holds across every level: its footprint, not its window.
u64 pyramid_entries(const plot::data::MinMaxPyramid& p) {
  u64 total = 0;
  for (u32 l = 0; l < p.level_count(); ++l) total += p.level_size(l);
  return total;
}

}  // namespace

PLOT_TEST(stream_series_sliding_window_queries_match_bruteforce) {
  // After the base has moved many times, an exact query must still match a scan
  // of exactly the retained samples. The pyramid's buckets live in absolute
  // sample space while the window's base is arbitrary, and that meeting point is
  // where an off-by-one hides -- so test the sliding case, not just the static
  // ArraySeries one.
  constexpr u64 kCapacity = 1u << 14;
  constexpr u64 kPushes = 60'000;
  const auto value = [](u64 i) { return f32((i * 29u) % 173u) - 86.0f; };
  plot::data::StreamSeries s(kCapacity, 0.0, 1e-6);
  for (u64 i = 0; i < kPushes; i += 337) {
    const u64 n = std::min<u64>(337, kPushes - i);
    for (u64 k = 0; k < n; ++k) s.push(value(i + k));
    s.acquire();   // the window's base moves on every one of these
  }

  const u64 n = s.size();
  const u64 base = s.produced() - n;
  PLOT_CHECK(n == kCapacity);
  PLOT_CHECK(base > 2 * kCapacity);             // the window really did slide
  PLOT_CHECK(s.pyramid().size() == n);          // pyramid() indexes the window
  PLOT_CHECK(s.pyramid().window_base() == base);
  PLOT_CHECK(s.pyramid().window_end() == s.produced());
  PLOT_CHECK_NEAR(s.x0(), f64(base) * 1e-6, 1e-12);
  PLOT_CHECK(s.x_bounds().min == s.x0());
  PLOT_CHECK_NEAR(s.x_bounds().max, f64(base + n - 1) * 1e-6, 1e-9);

  // Brute force in window coordinates over the retained samples only.
  const auto brute = [&](u64 a, u64 b) {
    f64 lo = f64(value(base + a)), hi = lo;
    f64 sum = 0.0, sq = 0.0;
    for (u64 i = a; i < b; ++i) {
      const f64 v = f64(value(base + i));
      lo = std::min(lo, v);
      hi = std::max(hi, v);
      sum += v;
      sq += v * v;
    }
    return std::array<f64, 4>{lo, hi, sum, sq};
  };

  const u64 starts[] = {0, 1, 63, 64, 65, 1000, 7919, n - 200};
  for (const u64 a : starts) {
    for (u64 b = a + 1; b <= n; b += 1301) {
      const std::array<f64, 4> want = brute(a, b);
      const plot::Range r = s.reduce(a, b);
      PLOT_CHECK_NEAR(r.min, want[0], 1e-4);
      PLOT_CHECK_NEAR(r.max, want[1], 1e-4);
      PLOT_CHECK_NEAR(s.sum(a, b), want[2], 1e-2);
      PLOT_CHECK_NEAR(s.sumsq(a, b), want[3], 1e-1);
    }
  }
  const std::array<f64, 4> whole = brute(0, n);
  PLOT_CHECK_NEAR(s.y_bounds().min, whole[0], 1e-6);
  PLOT_CHECK_NEAR(s.y_bounds().max, whole[1], 1e-6);

  // The level accessors are in absolute bucket coordinates: a stored bucket must
  // aggregate its own absolute span, the sliding window included.
  const u64 span = u64{kBucket} << 1;
  const u64 lbase = s.pyramid().level_base(1);
  const u64 lsize = s.pyramid().level_size(1);
  PLOT_CHECK(lsize >= 2);
  for (u64 b = lbase; b < lbase + std::min<u64>(lsize, 4); ++b) {
    const std::array<f64, 4> want = brute(b * span - base, (b + 1) * span - base);
    const plot::data::MinMax m = s.pyramid().bucket(1, b);
    PLOT_CHECK_NEAR(m.min, want[0], 1e-6);
    PLOT_CHECK_NEAR(m.max, want[1], 1e-6);
  }

  // level-1 gather: the envelope pairs of the retained buckets, bracketed by the
  // two partial edge buckets as raw samples, so the geometry spans the range
  // exactly. This is the mapping the renderer draws from. Bucket selection is
  // ceil/floor on the range edges (B3): `first` is the first level-1 bucket
  // wholly at or after window index 0, `last` one past the last wholly inside
  // the window -- so every emitted bucket's samples lie inside the window.
  std::vector<f32> out;
  s.gather(1, 0, n, out);
  const u64 first =
      std::max(base / span + (base % span == 0 ? u64{0} : u64{1}), lbase);   // ceil, then clamp
  const u64 last = std::min((base + n) / span, lbase + lsize);
  const u64 lead = first * span - base;              // raw pairs before the first whole bucket
  const u64 trail = last * span - base;              // first window index after the last one
  const u64 tail = trail < n ? n - trail : 0;
  const u64 env = last - first;
  PLOT_CHECK(out.size() == 2 * (lead + 2 * env + tail));   // each bucket is two pairs
  for (u64 i = 0; i < lead; ++i) {   // the partial leading bucket, raw
    // x relative to x0(): window index i emits i * dx, not the absolute
    // (base + i) * dx (relative-x contract).
    PLOT_CHECK_NEAR(out[2 * i], f32(f64(i) * 1e-6), 1e-9);
    PLOT_CHECK_NEAR(out[2 * i + 1], value(base + i), 1e-6);
  }
  for (u64 j = 0; j < env; ++j) {
    // Each whole bucket contributes two pairs: (first sample, min) and (last
    // sample, max). x is relative to x0(): window index i plots at i * dx.
    const u64 bi = first + j;
    const std::array<f64, 4> want = brute(bi * span - base, (bi + 1) * span - base);
    const u64 first_sample = bi * span - base;
    const u64 last_sample = std::min((bi + 1) * span, base + n) - 1 - base;
    const plot::usize lo = static_cast<plot::usize>(2 * (lead + 2 * j));
    const plot::usize hi = lo + 2;
    PLOT_CHECK_NEAR(out[lo], f32(f64(first_sample) * 1e-6), 1e-9);
    PLOT_CHECK_NEAR(out[lo + 1], want[0], 1e-6);
    PLOT_CHECK_NEAR(out[hi], f32(f64(last_sample) * 1e-6), 1e-9);
    PLOT_CHECK_NEAR(out[hi + 1], want[1], 1e-6);
  }
  for (u64 i = 0; i < tail; ++i) {   // the partial trailing bucket: the newest samples
    const plot::usize k = static_cast<plot::usize>(2 * (lead + 2 * env + i));
    PLOT_CHECK_NEAR(out[k], f32(f64(trail + i) * 1e-6), 1e-9);
    PLOT_CHECK_NEAR(out[k + 1], value(base + trail + i), 1e-6);
  }
}

PLOT_TEST(stream_series_steady_state_acquire_is_proportional_to_produced) {
  // The defect this pins: a saturated ring whose base moved every frame used to
  // rebuild the whole window, costing O(capacity * live_levels) merges per
  // acquire. It must cost O(samples produced since the last acquire) instead, so
  // the same production through a 16x larger window must cost the same.
  const auto merges_for = [](u64 capacity) {
    plot::data::StreamSeries s(capacity, 0.0, 1.0);
    constexpr u64 kSaturate = 1u << 20;
    for (u64 i = 0; i < kSaturate; ++i) s.push(f32(i % 977));
    s.acquire();
    PLOT_CHECK(s.size() == capacity);

    constexpr u64 kFrames = 64;
    constexpr u64 kFrame = 1u << 13;   // 8192 samples: 1 MS/s at ~122 fps
    plot::data::detail::merge_count().store(0, std::memory_order_relaxed);
    for (u64 f = 0; f < kFrames; ++f) {
      for (u64 k = 0; k < kFrame; ++k) s.push(f32((kSaturate + f * kFrame + k) % 977));
      s.acquire();
    }
    return plot::data::detail::merge_count().load(std::memory_order_relaxed);
  };

  constexpr u64 kProduced = 64 * (1u << 13);
  const u64 small = merges_for(1u << 16);   // 64k-sample window
  const u64 big = merges_for(1u << 20);     // 1M-sample window: 16x to rebuild
  // ~kProduced/32 in practice; rebuilding the 1M window alone would merge
  // kCapacity/64 * live_levels ~ 245'000 buckets per acquire.
  PLOT_CHECK(small <= kProduced / 16);
  PLOT_CHECK(big <= kProduced / 16);
  PLOT_CHECK(big <= 2 * small + kBucket);
}

PLOT_TEST(stream_series_pyramid_memory_is_bounded_by_the_window) {
  // 50M samples through a 1M-sample ring: the pyramid's footprint must track the
  // retained window (n/64 level-0 buckets plus the levels over them), not the
  // samples ever pushed.
  constexpr u64 kCapacity = 1u << 20;
  constexpr u64 kTotal = 50u * (1u << 20);
  constexpr u64 kChunk = 1u << 19;
  plot::data::StreamSeries s(kCapacity, 0.0, 1.0);
  u64 entries_after_10m = 0;
  for (u64 p = 0; p < kTotal; p += kChunk) {
    for (u64 i = 0; i < kChunk; ++i) s.push(f32((p + i) % 251));
    s.acquire();
    if (p + kChunk == 10u * (1u << 20)) entries_after_10m = pyramid_entries(s.pyramid());
  }
  const u64 entries = pyramid_entries(s.pyramid());
  const u64 bound = 2 * (kCapacity / kBucket) + kBucket;
  PLOT_CHECK(s.produced() == kTotal);
  PLOT_CHECK(s.size() == kCapacity);
  PLOT_CHECK(entries <= bound);
  // 5x the samples pushed, same footprint: levels deeper than one chunk can
  // differ by a bucket each, hence the slack.
  PLOT_CHECK(entries <= entries_after_10m + kBucket);
}

PLOT_TEST(gather_at_a_lod_level_reaches_both_ends_of_the_window) {
  // The pyramid stores complete buckets only, so its envelope stops at bucket
  // boundaries. Without the two partial edge buckets drawn from raw samples the
  // geometry would stop short of the newest samples -- the ones a live trace is
  // watched for -- and of the oldest retained ones.
  constexpr u64 kCapacity = 1u << 13;
  constexpr u64 kPushes = 40013;
  constexpr u64 kBase = kPushes - kCapacity;   // 31821: not a multiple of kBucket
  const auto value = [](u64 i) {
    if (i == kBase) return 25.0f;              // the oldest retained sample
    if (i == kPushes - 1) return -25.0f;       // the newest retained sample
    return f32(i % 97) * 0.25f - 12.0f;        // inside [-12, 12], so a bucket's
  };                                           // envelope can never equal either spike

  plot::data::StreamSeries s(kCapacity, 0.0, 1e-6);
  for (u64 i = 0; i < kPushes; ++i) s.push(value(i));
  s.acquire();
  const u64 n = s.size();
  PLOT_CHECK(n == kCapacity);
  PLOT_CHECK((kBase % kBucket) != 0);                       // both edges off the grid
  PLOT_CHECK((kPushes % (u64{kBucket} << 1)) != 0);

  const auto check_ends = [&](u32 level) {
    std::vector<f32> out;
    s.gather(level, 0, n, out);
    PLOT_CHECK(out.size() >= 8u);
    // Oldest: the first emitted pair is the retained sample at index 0, not the
    // first whole bucket (which starts `lead` samples later). x is relative to
    // x0(), so it emits 0 (relative-x contract).
    PLOT_CHECK_NEAR(out[0], 0.0, 1e-9);
    PLOT_CHECK_NEAR(out[1], 25.0, 1e-6);
    // Newest: the last emitted pair is the sample at index n-1, not the last
    // complete bucket.
    PLOT_CHECK_NEAR(out[out.size() - 2], f64(n - 1) * 1e-6, 1e-9);
    PLOT_CHECK_NEAR(out[out.size() - 1], -25.0, 1e-6);
  };
  PLOT_CHECK(s.pyramid().level_size(1) >= 2);
  check_ends(1);
  const u32 dense = s.pyramid().level_for(n, 1024.0);
  PLOT_CHECK(dense > 0u);
  check_ends(dense);

  // A static series whose window starts on the grid: only the trailing partial
  // bucket is at stake, and that is exactly what level_for() picks when a pixel
  // covers hundreds of samples.
  std::vector<f32> y(5000);
  for (u64 i = 0; i < y.size(); ++i) y[i] = value(i);
  y[0] = 25.0f;
  y[y.size() - 1] = -25.0f;
  const u64 ny = static_cast<u64>(y.size());
  PLOT_CHECK((ny % 256u) != 0u);   // 5000 = 19 * 256 + 136 at level 2
  plot::data::ArraySeries a(y, 0.0, 1e-6);
  std::vector<f32> out;
  a.gather(2, 0, ny, out);
  PLOT_CHECK(out.size() >= 8u);
  PLOT_CHECK_NEAR(out[0], 0.0, 1e-9);                       // level 2 starts on the grid
  PLOT_CHECK_NEAR(out[out.size() - 2], f32(f64(ny - 1) * 1e-6), 1e-9);
  PLOT_CHECK_NEAR(out[out.size() - 1], -25.0, 1e-6);
}

// --- B3: gather never emits geometry from outside [begin, end) ---------------

PLOT_TEST(neighbours_exact_hits_collapse_to_one_index) {
  // The contract says the bracket degenerates to (k, k) when x lands exactly on
  // a sample -- at the first sample too, where the old floor-only path answered
  // (0, 1) -- and stays (k, k+1) strictly between two samples.
  plot::data::ArraySeries s(std::vector<f32>{1, 2, 3, 4, 5}, 10.0, 0.25);
  u64 i0 = 9, i1 = 9;
  PLOT_CHECK(s.neighbours(s.x_at(0), i0, i1));    // first: was (0, 1)
  PLOT_CHECK(i0 == 0 && i1 == 0);
  PLOT_CHECK(s.neighbours(s.x_at(2), i0, i1));    // middle
  PLOT_CHECK(i0 == 2 && i1 == 2);
  PLOT_CHECK(s.neighbours(s.x_at(4), i0, i1));    // last
  PLOT_CHECK(i0 == 4 && i1 == 4);
  PLOT_CHECK(s.neighbours(s.x_at(1) + 0.1, i0, i1));   // strictly between
  PLOT_CHECK(i0 == 1 && i1 == 2);

  // A step whose k*dx does not round-trip through the division: the exact-hit
  // test must compare samples, not trust t's fraction being an integer.
  plot::data::ArraySeries t(std::vector<f32>{0, 1, 2, 3, 4, 5, 6, 7}, -1.0, 0.1);
  for (u64 i = 0; i < 8; ++i) {
    PLOT_CHECK(t.neighbours(t.x_at(i), i0, i1));
    PLOT_CHECK(i0 == i && i1 == i);
  }
  PLOT_CHECK(t.neighbours((t.x_at(3) + t.x_at(4)) / 2.0, i0, i1));
  PLOT_CHECK(i0 == 3 && i1 == 4);
}

PLOT_TEST(neighbours_outside_and_empty_are_false) {
  plot::data::ArraySeries s(std::vector<f32>{1, 2, 3}, 0.0, 1.0);
  u64 i0 = 7, i1 = 7;
  PLOT_CHECK(!s.neighbours(-0.5, i0, i1));    // left of the window
  PLOT_CHECK(!s.neighbours(3.5, i0, i1));     // right of the window
  PLOT_CHECK(!s.neighbours(std::nan(""), i0, i1));
  PLOT_CHECK(i0 == 7 && i1 == 7);             // untouched on false
  plot::data::ArraySeries e({}, 0.0, 1.0);
  PLOT_CHECK(!e.neighbours(1.0, i0, i1));
  PLOT_CHECK(e.snap_index(1.0) == 0);         // empty series snaps to 0
}

PLOT_TEST(gather_envelope_stays_inside_an_unaligned_range) {
  // The defect this pins: the first bucket was floor-selected and the bucket
  // containing `end - 1` was emitted whole, so an unaligned range drew
  // aggregates covering samples outside it. The samples outside [begin, end)
  // here are extreme on purpose: any envelope value beyond the exact reduce is
  // geometry that leaked in from outside the range.
  constexpr u64 kN = 4096;
  std::vector<f32> y(kN);
  for (u64 i = 0; i < kN; ++i) {
    y[i] = i < 1000 || i >= 3000 ? (i % 2 == 0 ? 1000.0f : -1000.0f)
                                 : 0.25f * f32(i % 4) - 0.375f;
  }
  plot::data::ArraySeries s(std::move(y), 0.0, 1.0);
  const u64 b = 1000 + 37, e = 3000 - 61;   // unaligned on every level's grid
  const plot::Range exact = s.reduce(b, e); // |values| <= ~0.375 inside the range
  PLOT_CHECK(exact.max < 1.0 && exact.min > -1.0);

  const auto check_range = [&](u32 level) {
    std::vector<f32> out;
    s.gather(level, b, e, out);
    PLOT_CHECK(out.size() >= 4u);
    for (plot::usize i = 0; i + 1 < out.size(); i += 2) {
      PLOT_CHECK(f64(out[i]) >= f64(b));    // emitted x inside [begin, end)
      PLOT_CHECK(f64(out[i]) < f64(e));
    }
    for (plot::usize i = 1; i < out.size(); i += 2) {
      PLOT_CHECK(f64(out[i]) >= exact.min - 1e-6);
      PLOT_CHECK(f64(out[i]) <= exact.max + 1e-6);
    }
  };
  check_range(1);
  check_range(2);
  check_range(3);
  check_range(6);

  // And the envelope values are not merely *inside* the range: each emitted
  // (min, max) pair equals the exact reduce over the samples that pair's two x
  // positions bracket. Walk the layout the contract fixes -- raw lead pairs,
  // two pairs per complete bucket, raw trail pairs -- so the raw edge spans
  // are not mistaken for envelope pairs.
  const u32 level = 2;
  const u64 span = u64{kBucket} << level;
  const u64 lend = s.pyramid().level_size(level);
  const u64 first = std::min(b / span + (b % span == 0 ? u64{0} : u64{1}), lend);
  const u64 last = std::min(e / span, lend);
  const u64 lead = first * span;
  const u64 buckets = last - first;
  PLOT_CHECK(lead > b && last * span < e);   // both edge spans really are partial
  std::vector<f32> out;
  s.gather(level, b, e, out);
  PLOT_CHECK(out.size() == 2 * ((lead - b) + 2 * buckets + (e - last * span)));
  plot::usize o = 0;
  for (u64 i = b; i < lead; ++i, o += 2) {   // the partial leading bucket, raw
    PLOT_CHECK_NEAR(out[o], f32(f64(i)), 1e-6);
    PLOT_CHECK(f64(out[o + 1]) >= exact.min - 1e-6 && f64(out[o + 1]) <= exact.max + 1e-6);
  }
  for (u64 j = 0; j < buckets; ++j) {
    const u64 bi = first + j;
    const u64 lo = bi * span;
    const u64 hi = std::min((bi + 1) * span, e) - 1;
    const plot::Range r = s.reduce(lo, hi + 1);
    PLOT_CHECK_NEAR(out[o], f32(f64(lo)), 1e-6);
    PLOT_CHECK_NEAR(out[o + 1], r.min, 1e-6);
    PLOT_CHECK_NEAR(out[o + 2], f32(f64(hi)), 1e-6);
    PLOT_CHECK_NEAR(out[o + 3], r.max, 1e-6);
    o += 4;
  }
  for (u64 i = last * span; i < e; ++i, o += 2) {   // the partial trailing bucket
    PLOT_CHECK_NEAR(out[o], f32(f64(i)), 1e-6);
    PLOT_CHECK(f64(out[o + 1]) >= exact.min - 1e-6 && f64(out[o + 1]) <= exact.max + 1e-6);
  }
  PLOT_CHECK(o == out.size());
}

PLOT_TEST(stream_gather_envelope_stays_inside_an_unaligned_range) {
  // Same invariant on a sliding window whose base is off the bucket grid. The
  // tame zone covers every sample the queried range reaches -- window [3808,
  // 12000) queried at window indices [300, n-77), absolute [4108, 11923) -- and
  // the extremes sit strictly before it, from 3808 to 4000.
  const auto value = [](u64 i) {
    return i < 4000 ? (i % 2 == 0 ? 500.0f : -500.0f)
                    : 0.125f * f32(i % 5) - 0.25f;
  };
  constexpr u64 kPushes = 12000;
  plot::data::StreamSeries s(1u << 13, 0.0, 1.0);   // window of 8192, base 3816
  for (u64 i = 0; i < kPushes; ++i) s.push(value(i));
  s.acquire();
  const u64 n = s.size();
  const u64 base = kPushes - n;
  PLOT_CHECK((base % kBucket) != 0);                // window starts mid-bucket
  const u64 b = 300, e = n - 77;                    // unaligned inside the window
  const plot::Range exact = s.reduce(b, e);
  PLOT_CHECK(exact.max < 1.0 && exact.min > -1.0);
  std::vector<f32> out;
  s.gather(1, b, e, out);
  PLOT_CHECK(out.size() >= 4u);
  for (plot::usize i = 0; i + 1 < out.size(); i += 2) {
    PLOT_CHECK(f64(out[i]) >= f64(b));              // relative x inside [begin, end)
    PLOT_CHECK(f64(out[i]) < f64(e));
  }
  for (plot::usize i = 1; i < out.size(); i += 2) {
    PLOT_CHECK(f64(out[i]) >= exact.min - 1e-6);
    PLOT_CHECK(f64(out[i]) <= exact.max + 1e-6);
  }
}

// --- B7: non-finite samples are excluded from every aggregate ----------------

PLOT_TEST(nan_samples_are_excluded_from_bounds_and_sums) {
  // NaN at the first, middle and last positions: bounds and reduce stay finite,
  // sums skip the NaN samples entirely.
  std::vector<f32> y{std::nanf(""), 1.0f, 2.0f, std::nanf(""), 3.0f, 4.0f,
                     std::nanf("")};
  plot::data::ArraySeries s(std::move(y), 0.0, 1.0);
  PLOT_CHECK_NEAR(s.y_bounds().min, 1.0, 1e-9);
  PLOT_CHECK_NEAR(s.y_bounds().max, 4.0, 1e-9);
  const plot::Range r = s.reduce(0, 7);
  PLOT_CHECK_NEAR(r.min, 1.0, 1e-9);
  PLOT_CHECK_NEAR(r.max, 4.0, 1e-9);
  PLOT_CHECK_NEAR(s.sum(0, 7), 10.0, 1e-9);
  PLOT_CHECK_NEAR(s.sumsq(0, 7), 30.0, 1e-9);
  // A range that is *all* NaN: zero-span, zero sums.
  PLOT_CHECK(s.reduce(0, 1).span() == 0.0);
  PLOT_CHECK(s.sum(0, 1) == 0.0);
  PLOT_CHECK(s.sumsq(0, 1) == 0.0);
  // And a NaN-only window reports a zero-span Range, not NaN.
  PLOT_CHECK(s.reduce(5, 7).span() == 0.0);
}

PLOT_TEST(stream_series_bounds_skip_nan_samples) {
  // acquire()'s bounds go through the same exact folds, so a stream with NaN
  // gaps publishes finite bounds -- and NaN at the edges of the window too.
  plot::data::StreamSeries s(256, 0.0, 1.0);
  s.push(std::nanf(""));                    // oldest retained sample
  for (int i = 0; i < 254; ++i) s.push(i % 2 == 0 ? -3.0f : 7.0f);
  s.push(std::nanf(""));                    // newest sample
  s.acquire();
  PLOT_CHECK(s.size() == 256);
  PLOT_CHECK_NEAR(s.y_bounds().min, -3.0, 1e-6);
  PLOT_CHECK_NEAR(s.y_bounds().max, 7.0, 1e-6);
  PLOT_CHECK_NEAR(s.reduce(0, s.size()).max, 7.0, 1e-6);
  PLOT_CHECK_NEAR(s.sum(0, s.size()), 127.0 * 4.0, 1e-3);
}

PLOT_TEST(window_without_finite_samples_aggregates_to_zero) {
  plot::data::ArraySeries s(std::vector<f32>(8, std::nanf("")), 0.0, 1.0);
  PLOT_CHECK(s.y_bounds().span() == 0.0);
  PLOT_CHECK(s.reduce(0, 8).span() == 0.0);
  PLOT_CHECK(s.sum(0, 8) == 0.0);
  PLOT_CHECK(s.sumsq(0, 8) == 0.0);
  // x bounds stay the grid's: the window exists even though nothing is finite.
  PLOT_CHECK_NEAR(s.x_bounds().min, 0.0, 1e-9);
  PLOT_CHECK_NEAR(s.x_bounds().max, 7.0, 1e-9);
}

// --- R3: any u32 gather level is safe ----------------------------------------

PLOT_TEST(gather_level_beyond_the_pyramid_clamps) {
  // R3: `kBucket << level` used to shift by >= 64 for a public level past 57
  // (undefined). Now any u32 level is accepted and clamped to the deepest one.
  std::vector<f32> y(8192);
  for (u64 i = 0; i < y.size(); ++i) y[i] = f32(i % 13) - 6.0f;
  plot::data::ArraySeries s(std::move(y), 0.0, 1.0);
  std::vector<f32> deepest;
  s.gather(plot::data::kMaxPyramidLevel, 0, s.size(), deepest);
  PLOT_CHECK(!deepest.empty());
  for (const u32 level : {u32{58}, u32{200}, u32{1000}, ~u32{0}}) {
    std::vector<f32> out;
    s.gather(level, 0, s.size(), out);
    PLOT_CHECK(out == deepest);   // clamped to the same deepest geometry
  }
  // The cap is where the span shift still fits: level 57 spans exactly 2^63.
  PLOT_CHECK(plot::data::kMaxPyramidLevel == 57);
  PLOT_CHECK((u64{kBucket} << plot::data::kMaxPyramidLevel) == (u64{1} << 63));
}

// --- R4: the ring's capacity is the rounded one -------------------------------

PLOT_TEST(stream_series_capacity_is_rounded_up_and_pre_sized) {
  // StreamSeries(100) retains 128 samples: the ring rounds up, size() saturates
  // at the rounded capacity, and the pyramid was reserved for it (a saturated
  // stream must not grow the level table mid-acquire).
  plot::data::StreamSeries s(100, 0.0, 1.0);
  for (int i = 0; i < 1000; ++i) s.push(f32(i));
  s.acquire();
  PLOT_CHECK(s.size() == 128);
  PLOT_CHECK_NEAR(s.y_bounds().max, 999.0, 1e-4);
  PLOT_CHECK_NEAR(s.x_bounds().max, s.x0() + 127.0, 1e-9);
  const plot::Range r = s.reduce(0, s.size());
  PLOT_CHECK_NEAR(r.min, 872.0, 1e-4);   // the oldest retained sample
}

// --- GatherMode: Envelope / Nearest / Mean ------------------------------------

PLOT_TEST(gather_modes_envelope_nearest_mean) {
  // 4096 samples of a sawtooth, level 2 (256-sample buckets, 16 of them).
  constexpr u64 kN = 4096;
  const auto value = [](u64 i) { return f32(i % 7) - 3.0f; };
  std::vector<f32> y(kN);
  for (u64 i = 0; i < kN; ++i) y[i] = value(i);
  plot::data::ArraySeries s(std::move(y), 0.0, 1.0);
  const u32 level = 2;
  const u64 span = u64{kBucket} << level;
  const u64 buckets = kN / span;
  PLOT_CHECK(buckets == 16);

  std::vector<f32> env, near, mean;
  s.gather(level, 0, kN, env, plot::data::GatherMode::Envelope);
  s.gather(level, 0, kN, near, plot::data::GatherMode::Nearest);
  s.gather(level, 0, kN, mean, plot::data::GatherMode::Mean);

  // Envelope: two points per bucket, min at the bucket's first sample, max at
  // its last (this sawtooth's per-bucket min is its first sample's value and
  // its max its last sample's value).
  PLOT_CHECK(env.size() == 4 * buckets);
  // Nearest and Mean: one point per bucket, x at the bucket's centre sample.
  PLOT_CHECK(near.size() == 2 * buckets);
  PLOT_CHECK(mean.size() == 2 * buckets);
  for (u64 j = 0; j < buckets; ++j) {
    const u64 lo = j * span;
    const u64 centre = lo + span / 2;
    const f64 centre_rel = f64(centre);   // dx == 1, relative to x0() == 0

    // x positions are relative and land where the mode says.
    PLOT_CHECK_NEAR(env[4 * j], f64(lo), 1e-6);
    PLOT_CHECK_NEAR(env[4 * j + 2], f64(lo + span - 1), 1e-6);
    PLOT_CHECK_NEAR(near[2 * j], centre_rel, 1e-6);
    PLOT_CHECK_NEAR(mean[2 * j], centre_rel, 1e-6);

    // Envelope y is the exact per-bucket min/max.
    const plot::Range r = s.reduce(lo, lo + span);
    PLOT_CHECK_NEAR(env[4 * j + 1], r.min, 1e-6);
    PLOT_CHECK_NEAR(env[4 * j + 3], r.max, 1e-6);
    // Nearest y is the raw centre sample.
    PLOT_CHECK_NEAR(near[2 * j + 1], value(centre), 1e-6);
    // Mean y is the bucket's mean.
    const f64 want_mean = s.sum(lo, lo + span) / f64(span);
    PLOT_CHECK_NEAR(mean[2 * j + 1], want_mean, 1e-6);
  }

  // Level 0 and sub-bucket ranges are the raw samples for every mode.
  std::vector<f32> raw;
  s.gather(0, 5, 300, raw, plot::data::GatherMode::Mean);
  PLOT_CHECK(raw.size() == 2 * (300 - 5));
  PLOT_CHECK_NEAR(raw[0], 5.0, 1e-9);
}

PLOT_TEST(gather_skips_all_nan_buckets) {
  // Level-1 bucket 1 (samples [128, 256)) is entirely NaN: the exact queries
  // treat it as absent, the gathered geometry emits nothing for it in Envelope
  // and Mean, and Nearest draws its centre sample raw -- NaN, as level 0 would.
  std::vector<f32> y(384, 1.0f);
  for (u64 i = 128; i < 256; ++i) y[i] = std::nanf("");
  plot::data::ArraySeries s(std::move(y), 0.0, 1.0);
  PLOT_CHECK(s.reduce(128, 256).span() == 0.0);
  PLOT_CHECK(s.sum(128, 256) == 0.0);
  const plot::Range r = s.reduce(0, 384);
  PLOT_CHECK_NEAR(r.min, 1.0, 1e-9);
  PLOT_CHECK_NEAR(r.max, 1.0, 1e-9);

  // first = 0, last = 256/128 = 2: buckets 0 and 1, no lead/trail raw spans.
  std::vector<f32> env, mean, near;
  s.gather(1, 0, 256, env, plot::data::GatherMode::Envelope);
  PLOT_CHECK(env.size() == 4);              // bucket 0's two pairs only
  PLOT_CHECK_NEAR(env[0], 0.0, 1e-9);
  PLOT_CHECK_NEAR(env[1], 1.0, 1e-9);
  PLOT_CHECK_NEAR(env[2], 127.0, 1e-9);
  PLOT_CHECK_NEAR(env[3], 1.0, 1e-9);
  s.gather(1, 0, 256, mean, plot::data::GatherMode::Mean);
  PLOT_CHECK(mean.size() == 2);             // bucket 1 has no finite mean: skipped
  PLOT_CHECK_NEAR(mean[0], 64.0, 1e-9);
  PLOT_CHECK_NEAR(mean[1], 1.0, 1e-9);
  s.gather(1, 0, 256, near, plot::data::GatherMode::Nearest);
  PLOT_CHECK(near.size() == 4);             // every bucket still has a centre sample
  PLOT_CHECK_NEAR(near[0], 64.0, 1e-9);
  PLOT_CHECK_NEAR(near[1], 1.0, 1e-9);
  PLOT_CHECK_NEAR(near[2], 192.0, 1e-9);    // bucket 1's centre, drawn raw
  PLOT_CHECK(std::isnan(near[3]));
}

// --- Snapshot: StreamSeries::freeze() -----------------------------------------

PLOT_TEST(stream_series_freeze_copies_the_last_acquire_window) {
  plot::data::StreamSeries s(1024, 5.0, 0.5);
  for (u64 i = 0; i < 400; ++i) s.push(f32(i % 31));
  s.acquire();
  const u64 n = s.size();
  const plot::Range bounds = s.y_bounds();
  std::vector<f32> gathered;
  s.gather(1, 0, n, gathered);
  const plot::Range reduced = s.reduce(0, n);

  // Freeze now, while the ring still holds exactly the published window.
  const plot::data::ArraySeries frozen = s.freeze();
  PLOT_CHECK(frozen.size() == n);
  PLOT_CHECK_NEAR(frozen.x0(), 5.0, 1e-12);
  PLOT_CHECK_NEAR(frozen.dx(), 0.5, 1e-12);
  PLOT_CHECK_NEAR(frozen.y_bounds().min, bounds.min, 1e-6);
  PLOT_CHECK_NEAR(frozen.y_bounds().max, bounds.max, 1e-6);
  const plot::Range fr = frozen.reduce(0, n);
  PLOT_CHECK_NEAR(fr.min, reduced.min, 1e-6);
  PLOT_CHECK_NEAR(fr.max, reduced.max, 1e-6);
  PLOT_CHECK_NEAR(frozen.sum(0, n), s.sum(0, n), 1e-3);
  std::vector<f32> frozen_gathered;
  frozen.gather(1, 0, n, frozen_gathered);
  PLOT_CHECK(frozen_gathered == gathered);   // same geometry, byte for byte

  // The producer now laps the ring (700 more pushes overwrite the oldest
  // slots) and a later acquire() grows the published window -- the frozen copy
  // answers from its own samples and never moves.
  for (u64 i = 0; i < 700; ++i) s.push(-1000.0f);
  s.acquire();
  PLOT_CHECK(s.size() == 1024);
  PLOT_CHECK(frozen.size() == n);
  PLOT_CHECK_NEAR(frozen.y_bounds().max, bounds.max, 1e-6);
  PLOT_CHECK_NEAR(frozen.y_bounds().min, bounds.min, 1e-6);
  PLOT_CHECK(frozen.reduce(0, n).max == reduced.max);
}

// --- Incremental append: ArraySeries::push() ----------------------------------

PLOT_TEST(array_series_push_matches_commit) {
  // The incremental path must land exactly where the wholesale path does:
  // bounds, x bounds, exact queries and the pyramid's buckets themselves.
  std::vector<f32> y;
  y.reserve(1000);
  for (u64 i = 0; i < 1000; ++i) {
    y.push_back(i % 97 == 0 ? std::nanf("") : f32((i * 37) % 251) - 125.0f);
  }
  plot::data::ArraySeries pushed({}, 10.0, 0.25);
  for (const f32 v : y) pushed.push(v);
  const plot::data::ArraySeries committed(y, 10.0, 0.25);

  PLOT_CHECK(pushed.size() == committed.size());
  PLOT_CHECK_NEAR(pushed.x_bounds().min, committed.x_bounds().min, 1e-12);
  PLOT_CHECK_NEAR(pushed.x_bounds().max, committed.x_bounds().max, 1e-12);
  PLOT_CHECK_NEAR(pushed.y_bounds().min, committed.y_bounds().min, 1e-6);
  PLOT_CHECK_NEAR(pushed.y_bounds().max, committed.y_bounds().max, 1e-6);
  for (u64 a = 0; a < 1000; a += 173) {
    for (u64 b = a + 1; b <= 1000; b += 331) {
      const plot::Range ra = pushed.reduce(a, b);
      const plot::Range rb = committed.reduce(a, b);
      PLOT_CHECK_NEAR(ra.min, rb.min, 1e-9);
      PLOT_CHECK_NEAR(ra.max, rb.max, 1e-9);
      PLOT_CHECK_NEAR(pushed.sum(a, b), committed.sum(a, b), 1e-9);
      PLOT_CHECK_NEAR(pushed.sumsq(a, b), committed.sumsq(a, b), 1e-9);
    }
  }
  for (u32 l = 0; l < 3; ++l) {
    PLOT_CHECK(pushed.pyramid().level_base(l) == committed.pyramid().level_base(l));
    PLOT_CHECK(pushed.pyramid().level_size(l) == committed.pyramid().level_size(l));
    for (u64 b = pushed.pyramid().level_base(l);
         b < pushed.pyramid().level_base(l) + pushed.pyramid().level_size(l); ++b) {
      const plot::data::MinMax ma = pushed.pyramid().bucket(l, b);
      const plot::data::MinMax mb = committed.pyramid().bucket(l, b);
      PLOT_CHECK_NEAR(ma.min, mb.min, 1e-9);
      PLOT_CHECK_NEAR(ma.max, mb.max, 1e-9);
      PLOT_CHECK_NEAR(pushed.pyramid().bucket_sum(l, b), committed.pyramid().bucket_sum(l, b),
                      1e-9);
      PLOT_CHECK(pushed.pyramid().bucket_count(l, b) == committed.pyramid().bucket_count(l, b));
    }
  }
}