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

#include <plot/data/xy_series.hpp>

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <utility>

namespace plot::data {
namespace {

/// Ceiling division: the first level-`l` bucket whose whole span starts at or
/// after `v` samples (in bucket units).
constexpr u64 ceil_div(u64 v, u64 d) { return v / d + (v % d == 0 ? u64{0} : u64{1}); }

}  // namespace

XYSeries::XYSeries(std::vector<f64> x, std::vector<f32> y) { reset(std::move(x), std::move(y)); }

void XYSeries::reset(std::vector<f64> x, std::vector<f32> y) {
  // The binary searches below (and every O(log n) consumer of them) are only
  // correct on sorted input; see the header for the precondition.
  assert(x.size() == y.size());
  assert(std::is_sorted(x.begin(), x.end()) && "XYSeries requires x sorted ascending");
  x_ = std::move(x);
  y_ = std::move(y);
  commit();
}

void XYSeries::commit() {
  const u64 n = static_cast<u64>(y_.size());
  pyramid_.clear();
  pyramid_.reserve(n);

  // Same NaN-aware fold as ArraySeries::commit(): non-finite y is excluded
  // from the cached bounds, and a window with no finite sample aggregates to a
  // zero-span Range. (The pyramid's append() applies the same policy per
  // bucket.) The x bounds come straight from the sorted array's ends.
  f32 lo = 0.0f, hi = 0.0f;
  bool finite = false;
  for (const f32 v : y_) {
    pyramid_.append(v);
    if (!std::isfinite(v)) continue;
    if (finite) {
      lo = std::min(lo, v);
      hi = std::max(hi, v);
    } else {
      lo = v;
      hi = v;
      finite = true;
    }
  }
  y_has_finite_ = finite;
  y_bounds_ = finite ? Range{f64(lo), f64(hi)} : Range{0.0, 0.0};
  x_bounds_ = n == 0 ? Range{0.0, 0.0} : Range{x_.front(), x_.back()};
}

u64 XYSeries::size() const { return static_cast<u64>(y_.size()); }

f64 XYSeries::x0() const { return x_.empty() ? 0.0 : x_.front(); }

f64 XYSeries::dx() const { return 0.0; }

Range XYSeries::x_bounds() const { return x_bounds_; }

Range XYSeries::y_bounds() const { return y_bounds_; }

const MinMaxPyramid& XYSeries::pyramid() const { return pyramid_; }

f64 XYSeries::x_at(u64 i) const { return x_[static_cast<usize>(i)]; }

bool XYSeries::neighbours(f64 x, u64& i0, u64& i1) const {
  const u64 n = size();
  if (n == 0) return false;
  // Also rejects NaN: every comparison with it is false.
  if (!(x >= x_.front() && x <= x_.back())) return false;
  // Last sample at or left of x, first at or right of it. Both searches are in
  // range: x >= front puts a lower_bound hit at index >= 1's predecessor at 0,
  // and x <= back caps upper_bound's predecessor at n - 1.
  const auto up = std::upper_bound(x_.begin(), x_.end(), x);
  const auto lo = std::lower_bound(x_.begin(), x_.end(), x);
  i0 = static_cast<u64>(up - x_.begin()) - 1;
  i1 = static_cast<u64>(lo - x_.begin());
  return true;
}

u64 XYSeries::snap_index(f64 x) const {
  const u64 n = size();
  if (n == 0) return 0;
  if (!(x > x_.front())) return 0;      // at or left of the first sample; NaN too
  if (x >= x_.back()) return n - 1;
  const auto it = std::lower_bound(x_.begin(), x_.end(), x);
  const u64 hi = static_cast<u64>(it - x_.begin());
  const u64 lo = hi - 1;                // hi >= 1: x is strictly past x_.front()
  return (x - x_[static_cast<usize>(lo)]) <= (x_[static_cast<usize>(hi)] - x)
             ? lo
             : hi;                      // ties go to the lower index
}

void XYSeries::gather(u32 level, u64 begin, u64 end, std::vector<f32>& out_xy,
                      GatherMode mode) const {
  const u64 n = size();
  const u64 b = std::min(begin, n);
  const u64 e = std::min(std::max(end, b), n);
  if (b >= e) return;

  // x emitted relative to x0(): see Series::gather(). The x array is f64, so
  // the subtraction happens wide and only the small, exact offset narrows.
  const f64 origin = x0();

  if (level == 0) {
    for (u64 i = b; i < e; ++i) {
      out_xy.push_back(f32(x_[static_cast<usize>(i)] - origin));
      out_xy.push_back(y_[static_cast<usize>(i)]);
    }
    return;
  }
  level = std::min(level, kMaxPyramidLevel);

  // The bucket grid is the sample index, exactly as for a uniform series, and
  // an XYSeries never slides its window, so the pyramid's base is 0 and the
  // window range needs no origin shift. Bucket selection is the B3-correct
  // ceil/floor: `first` is the first bucket whose samples all sit at or after
  // `begin`, `last` one past the last bucket wholly before `end`, and the two
  // edge spans around them are read from the raw arrays -- so no emitted
  // geometry comes from outside [begin, end), and the emitted x values are the
  // array's own coordinates relative to x0().
  const u64 span = u64{kBucket} << level;
  const u64 lend = pyramid().level_size(level);
  const u64 first = std::min(ceil_div(b, span), lend);
  const u64 last = std::min(e / span, lend);
  if (first >= last) {
    for (u64 i = b; i < e; ++i) {
      out_xy.push_back(f32(x_[static_cast<usize>(i)] - origin));
      out_xy.push_back(y_[static_cast<usize>(i)]);
    }
    return;
  }
  const u64 lead = first * span;
  for (u64 i = b; i < lead; ++i) {
    out_xy.push_back(f32(x_[static_cast<usize>(i)] - origin));
    out_xy.push_back(y_[static_cast<usize>(i)]);
  }

  switch (mode) {
    case GatherMode::Envelope: {
      for (u64 bi = first; bi < last; ++bi) {
        if (pyramid().bucket_count(level, bi) == 0) continue;
        const MinMax m = pyramid().bucket(level, bi);
        const u64 lo = bi * span;                          // the bucket's first sample
        const u64 hi = std::min((bi + 1) * span, e) - 1;   // its last sample inside the range
        out_xy.push_back(f32(x_[static_cast<usize>(lo)] - origin));
        out_xy.push_back(m.min);
        out_xy.push_back(f32(x_[static_cast<usize>(hi)] - origin));
        out_xy.push_back(m.max);
      }
      break;
    }
    case GatherMode::Nearest: {
      for (u64 bi = first; bi < last; ++bi) {
        const u64 centre = bi * span + span / 2;   // inside the bucket, hence inside [b, e)
        out_xy.push_back(f32(x_[static_cast<usize>(centre)] - origin));
        out_xy.push_back(y_[static_cast<usize>(centre)]);
      }
      break;
    }
    case GatherMode::Mean: {
      for (u64 bi = first; bi < last; ++bi) {
        const u32 count = pyramid().bucket_count(level, bi);
        if (count == 0) continue;
        const u64 centre = bi * span + span / 2;
        out_xy.push_back(f32(x_[static_cast<usize>(centre)] - origin));
        out_xy.push_back(f32(pyramid().bucket_sum(level, bi) / f64(count)));
      }
      break;
    }
  }

  const u64 trail = last * span;
  for (u64 i = trail; i < e; ++i) {
    out_xy.push_back(f32(x_[static_cast<usize>(i)] - origin));
    out_xy.push_back(y_[static_cast<usize>(i)]);
  }
}

void XYSeries::copy_raw(u64 begin, u64 end, f32* out) const {
  const u64 want = end > begin ? end - begin : 0;
  if (want == 0) return;
  const u64 n = size();
  const u64 b = std::min(begin, n);
  const u64 e = std::min(std::max(end, b), n);
  std::copy(y_.begin() + static_cast<std::ptrdiff_t>(b),
            y_.begin() + static_cast<std::ptrdiff_t>(e), out);
  std::fill(out + static_cast<std::ptrdiff_t>(e - b),
            out + static_cast<std::ptrdiff_t>(want), 0.0f);
}

}  // namespace plot::data