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

#include "harness.hpp"
#include <plot/data/array_series.hpp>
#include <plot/data/xy_series.hpp>
#include <algorithm>
#include <cmath>
#include <random>
#include <vector>

using plot::data::kBucket;

namespace {

/// Monotone but decidedly non-uniform x: early samples a millisecond apart,
/// later ones seconds apart. Exact in f64, no fixed dx anywhere.
std::vector<f64> spread_x(u64 n) {
  std::vector<f64> x;
  x.reserve(static_cast<plot::usize>(n));
  f64 v = 0.0;
  for (u64 i = 0; i < n; ++i) {
    x.push_back(v);
    v += 0.001 + f64(i) * 0.01;
  }
  return x;
}

std::vector<f32> saw_y(u64 n) {
  std::vector<f32> y;
  y.reserve(static_cast<plot::usize>(n));
  for (u64 i = 0; i < n; ++i) y.push_back(f32((i * 13) % 97) - 48.0f);
  return y;
}

}  // namespace

PLOT_TEST(xy_series_reports_its_own_grid) {
  auto x = spread_x(500);
  const std::vector<f64> xs = x;   // keep an un-moved copy for the checks
  plot::data::XYSeries s(std::move(x), saw_y(500));
  PLOT_CHECK(s.size() == 500);
  PLOT_CHECK_NEAR(s.x0(), xs.front(), 1e-12);
  PLOT_CHECK(s.dx() == 0.0);       // meaningless by contract; nothing reads it
  for (const u64 i : {u64{0}, u64{1}, u64{250}, u64{499}}) {
    PLOT_CHECK_NEAR(s.x_at(i), xs[static_cast<plot::usize>(i)], 1e-12);
  }
  // x_bounds come from the array's ends, wherever they sit on the number line.
  PLOT_CHECK_NEAR(s.x_bounds().min, xs.front(), 1e-12);
  PLOT_CHECK_NEAR(s.x_bounds().max, xs.back(), 1e-12);
  PLOT_CHECK_NEAR(s.y_bounds().min, -48.0, 1e-6);
  PLOT_CHECK_NEAR(s.y_bounds().max, 48.0, 1e-6);

  // Epoch-style offset: the same shape shifted by 1.7e9 keeps exact bounds, and
  // gather() emits x relative to x0() -- small numbers, not quantised epochs.
  std::vector<f64> epoch = spread_x(500);
  const std::vector<f64> es = epoch;   // un-moved copy for the expectations
  for (auto& v : epoch) v += 1.7e9;
  plot::data::XYSeries e(std::move(epoch), saw_y(500));
  PLOT_CHECK_NEAR(e.x_bounds().min, 1.7e9, 1e-3);
  PLOT_CHECK_NEAR(e.x_bounds().max, 1.7e9 + (xs.back() - xs.front()), 1e-3);
  std::vector<f32> out;
  e.gather(0, 0, 500, out);
  PLOT_CHECK(out.size() == 1000);
  // x is relative to x0(): x_at(i) - x0(), the array's own offsets -- small
  // values even though the absolute epoch is 1.7e9. The expectations come from
  // the epoch array itself, since that is the subtraction the contract does.
  PLOT_CHECK_NEAR(out[0], 0.0, 1e-9);                          // x_at(0) - x0()
  PLOT_CHECK_NEAR(out[2], f32(es[1] - es[0]), 1e-6);
  PLOT_CHECK_NEAR(out[998], f32(es[499] - es[0]), 1e-3);       // span, not 1.7e9
  PLOT_CHECK(out[998] < 2000.0f);
}

PLOT_TEST(xy_series_neighbour_lookups_crosshair_style) {
  // The crosshair asks neighbours() at an arbitrary x and reports whichever of
  // the two samples it gets back is nearer. On a non-uniform grid that needs
  // real searches: a uniform-grid formula would land megabytes off here.
  auto x = spread_x(1000);
  const std::vector<f64> xs = x;
  plot::data::XYSeries s(std::move(x), saw_y(1000));
  const auto at = [&](u64 i) { return xs[static_cast<plot::usize>(i)]; };

  u64 i0 = 7, i1 = 7;
  // Midway between samples 100 and 101: brackets those two.
  const f64 mid = (at(100) + at(101)) / 2.0;
  PLOT_CHECK(s.neighbours(mid, i0, i1));
  PLOT_CHECK(i0 == 100);
  PLOT_CHECK(i1 == 101);
  // Exactly on a sample: both indices are that sample.
  PLOT_CHECK(s.neighbours(at(250), i0, i1));
  PLOT_CHECK(i0 == 250 && i1 == 250);
  // Just left of a sample: still bracketed from below.
  PLOT_CHECK(s.neighbours(std::nextafter(at(250), 0.0), i0, i1));
  PLOT_CHECK(i0 == 249 && i1 == 250);
  // Edges are inclusive.
  PLOT_CHECK(s.neighbours(at(0), i0, i1));
  PLOT_CHECK(i0 == 0 && i1 == 0);
  PLOT_CHECK(s.neighbours(at(999), i0, i1));
  PLOT_CHECK(i0 == 999 && i1 == 999);   // exact hit on the last sample
  // Outside the bounds (and NaN): false, indices untouched.
  const u64 keep = 42;
  i0 = i1 = keep;
  PLOT_CHECK(!s.neighbours(at(0) - 1.0, i0, i1));
  PLOT_CHECK(!s.neighbours(at(999) + 1.0, i0, i1));
  PLOT_CHECK(!s.neighbours(std::nan(""), i0, i1));
  PLOT_CHECK(i0 == keep && i1 == keep);
  // Clustered duplicates: the contract is "last index <= x" and "first index
  // >= x", so a run of equal x brackets the run's ends.
  std::vector<f64> dup{1.0, 2.0, 2.0, 2.0, 3.0};
  plot::data::XYSeries d(dup, {1.0f, 2.0f, 3.0f, 4.0f, 5.0f});
  PLOT_CHECK(d.neighbours(2.0, i0, i1));
  PLOT_CHECK(i0 == 3 && i1 == 1);   // last <= 2 is index 3, first >= 2 is index 1
}

PLOT_TEST(xy_series_snap_index_finds_the_nearest_sample) {
  auto x = spread_x(1000);
  const std::vector<f64> xs = x;
  plot::data::XYSeries s(std::move(x), saw_y(1000));
  const auto at = [&](u64 i) { return xs[static_cast<plot::usize>(i)]; };

  PLOT_CHECK(s.snap_index(at(400)) == 400);                       // exact hit
  PLOT_CHECK(s.snap_index((at(400) + at(401)) / 2.0 + 1e-9) == 401);
  PLOT_CHECK(s.snap_index((at(400) + at(401)) / 2.0 - 1e-9) == 400);
  PLOT_CHECK(s.snap_index(at(0) - 100.0) == 0);                   // clamped
  PLOT_CHECK(s.snap_index(at(999) + 100.0) == 999);               // clamped
  PLOT_CHECK(s.snap_index(std::nan("")) == 0);                    // NaN lands on 0
  // A tie goes to the lower index (midpoints are impossible on this x, so make
  // one explicitly).
  std::vector<f64> tie{0.0, 10.0};
  plot::data::XYSeries t(tie, {1.0f, 2.0f});
  PLOT_CHECK(t.snap_index(5.0) == 0);
  PLOT_CHECK(t.snap_index(5.0 + 1e-9) == 1);
}

PLOT_TEST(xy_series_pyramid_matches_array_series_over_the_same_y) {
  // The pyramid buckets by sample index, exactly as for a uniform series, so an
  // XYSeries over x = i*dx must answer every exact query identically to the
  // ArraySeries over the same y -- and the bucket aggregates must agree too.
  constexpr u64 kN = 3000;
  std::mt19937 rng(11);
  std::uniform_real_distribution<float> d(-50.0f, 50.0f);
  std::vector<f32> y(kN);
  for (auto& v : y) v = d(rng);
  std::vector<f64> x(kN);
  for (u64 i = 0; i < kN; ++i) x[i] = f64(i) * 0.5;
  const plot::data::XYSeries xy(x, y);
  const plot::data::ArraySeries uni(y, 0.0, 0.5);

  PLOT_CHECK_NEAR(xy.x0(), uni.x0(), 1e-12);
  for (u64 a = 0; a < kN; a += 211) {
    for (u64 b = a + 1; b <= kN; b += 499) {
      const plot::Range rx = xy.reduce(a, b);
      const plot::Range ru = uni.reduce(a, b);
      PLOT_CHECK_NEAR(rx.min, ru.min, 1e-6);
      PLOT_CHECK_NEAR(rx.max, ru.max, 1e-6);
      PLOT_CHECK_NEAR(xy.sum(a, b), uni.sum(a, b), 1e-6);
      PLOT_CHECK_NEAR(xy.sumsq(a, b), uni.sumsq(a, b), 1e-6);
    }
  }
  for (u32 l = 0; l < 3; ++l) {
    PLOT_CHECK(xy.pyramid().level_base(l) == uni.pyramid().level_base(l));
    PLOT_CHECK(xy.pyramid().level_size(l) == uni.pyramid().level_size(l));
    for (u64 b = xy.pyramid().level_base(l);
         b < xy.pyramid().level_base(l) + xy.pyramid().level_size(l); ++b) {
      PLOT_CHECK_NEAR(xy.pyramid().bucket(l, b).min, uni.pyramid().bucket(l, b).min, 1e-6);
      PLOT_CHECK_NEAR(xy.pyramid().bucket(l, b).max, uni.pyramid().bucket(l, b).max, 1e-6);
    }
  }
}

PLOT_TEST(xy_series_gather_envelope_matches_exact_reduces) {
  // gather() buckets by index with the same B3 ceil/floor selection, but the x
  // of every emitted point is the array's own coordinate relative to x0() --
  // and each envelope pair equals the exact reduce over the samples its two x
  // positions bracket.
  constexpr u64 kN = 4096;
  auto x = spread_x(kN);
  const std::vector<f64> xs = x;
  const std::vector<f32> ys = saw_y(kN);
  plot::data::XYSeries s(std::move(x), ys);
  const auto at = [&](u64 i) { return xs[static_cast<plot::usize>(i)]; };
  const auto val = [&](u64 i) { return ys[static_cast<plot::usize>(i)]; };
  const f64 x0 = at(0);

  const u32 level = 2;
  const u64 span = u64{kBucket} << level;
  const u64 b = 500 + 33, e = 3500 - 57;   // unaligned on the bucket grid
  std::vector<f32> out;
  s.gather(level, b, e, out, plot::data::GatherMode::Envelope);
  PLOT_CHECK(out.size() >= 8u);
  // Emitted x stays inside [x_at(begin), x_at(end)) and the geometry reaches
  // both the range's extremes (the lead/trail raw spans carry them).
  const plot::Range exact = s.reduce(b, e);
  f64 emitted_min = 1e30, emitted_max = -1e30;
  for (plot::usize i = 0; i + 1 < out.size(); i += 2) {
    const f64 rel = f64(out[i]);
    // The emitted x is f32(x_at(i) - x0()), so allow the f32 cast's slack
    // (relative positions here reach ~1.4e4, where f32 ULP is ~1e-3).
    PLOT_CHECK(rel >= at(b) - x0 - 1e-3);
    PLOT_CHECK(rel <= at(e - 1) - x0 + 1e-3);
    emitted_min = std::min(emitted_min, f64(out[i + 1]));
    emitted_max = std::max(emitted_max, f64(out[i + 1]));
  }
  PLOT_CHECK_NEAR(emitted_min, exact.min, 1e-6);
  PLOT_CHECK_NEAR(emitted_max, exact.max, 1e-6);
  // Envelope pairs == exact reduce over their own spans. Walk the layout the
  // gather contract fixes -- raw lead pairs, then two pairs per complete
  // bucket, then raw trail pairs -- and read each pair's span off its x
  // positions (all saw_y values are finite, so no bucket is skipped).
  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(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(at(i) - x0), 1e-3);
    PLOT_CHECK_NEAR(out[o + 1], val(i), 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(at(lo) - x0), 1e-3);
    PLOT_CHECK_NEAR(out[o + 1], r.min, 1e-6);
    PLOT_CHECK_NEAR(out[o + 2], f32(at(hi) - x0), 1e-3);
    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(at(i) - x0), 1e-3);
    PLOT_CHECK_NEAR(out[o + 1], val(i), 1e-6);
  }
  PLOT_CHECK(o == out.size());
}

PLOT_TEST(xy_series_nan_policy_matches_the_other_series) {
  // Non-finite y: excluded from bounds and sums; an all-NaN window aggregates
  // to zero-span / zero.
  std::vector<f64> x{0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
  plot::data::XYSeries s(x, {std::nanf(""), 1.0f, 2.0f, std::nanf(""), 3.0f, std::nanf("")});
  PLOT_CHECK_NEAR(s.y_bounds().min, 1.0, 1e-9);
  PLOT_CHECK_NEAR(s.y_bounds().max, 3.0, 1e-9);
  const plot::Range r = s.reduce(0, 6);
  PLOT_CHECK_NEAR(r.min, 1.0, 1e-9);
  PLOT_CHECK_NEAR(r.max, 3.0, 1e-9);
  PLOT_CHECK_NEAR(s.sum(0, 6), 6.0, 1e-9);
  PLOT_CHECK(s.reduce(4, 5).span() == 0.0);
  plot::data::XYSeries all(x, std::vector<f32>(6, std::nanf("")));
  PLOT_CHECK(all.y_bounds().span() == 0.0);
  PLOT_CHECK(all.sum(0, 6) == 0.0);
  PLOT_CHECK_NEAR(all.x_bounds().min, 0.0, 1e-12);   // x bounds are the array's
  PLOT_CHECK_NEAR(all.x_bounds().max, 5.0, 1e-12);
}

PLOT_TEST(xy_series_commit_rebuilds_and_reset_replaces) {
  auto x = spread_x(256);
  plot::data::XYSeries s(std::move(x), saw_y(256));
  const plot::Range before = s.y_bounds();
  std::vector<f64> x2{0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0, 32.0};
  s.reset(std::move(x2), {5.0f, 3.0f, 1.0f, 4.0f, 4.0f, 4.0f, 4.0f, 4.0f});
  PLOT_CHECK(s.size() == 8);
  PLOT_CHECK_NEAR(s.y_bounds().min, 1.0, 1e-9);
  PLOT_CHECK_NEAR(s.y_bounds().max, 5.0, 1e-9);
  PLOT_CHECK_NEAR(s.x_bounds().max, 32.0, 1e-12);
  PLOT_CHECK(s.y_bounds().max != before.max);
  // reset() re-asserts the precondition too: the arrays above are sorted, so
  // this passes; unsorted x is a debug-abort by contract (see the header) and
  // is deliberately not exercised here.
}