// Numerically stable reductions (FP64 reference implementations).
// sum_naive is provided only to *measure* error; library code uses sum_neumaier or sum_pairwise.
#pragma once
#include <span>
#include <cmath>
#include <vector>
#include <algorithm>
#include "raid/core/types.hpp"
#include "raid/core/error.hpp"

namespace raid::stats {

template <class T>
T sum_naive(std::span<const T> x) { T s = 0; for (T v : x) s += v; return s; }

/// Neumaier's improved Kahan–Babuška summation.
template <class T>
T sum_neumaier(std::span<const T> x) {
  T s = 0, c = 0;
  for (T v : x) {
    const T t = s + v;
    if (std::abs(s) >= std::abs(v)) c += (s - t) + v; else c += (v - t) + s;
    s = t;
  }
  return s + c;
}

/// Pairwise (cascade) summation with a base-case block of 128; error O(log n) ulp.
template <class T>
T sum_pairwise(std::span<const T> x) {
  if (x.size() <= 128) { T s = 0; for (T v : x) s += v; return s; }
  const usize m = x.size() / 2;
  return sum_pairwise(x.subspan(0, m)) + sum_pairwise(x.subspan(m));
}

template <class T>
T dot_neumaier(std::span<const T> a, std::span<const T> b) {
  RAID_REQUIRE(a.size() == b.size(), invalid_argument, "dot: size mismatch");
  T s = 0, c = 0;
  for (usize i = 0; i < a.size(); ++i) {
    const T v = a[i] * b[i], t = s + v;
    if (std::abs(s) >= std::abs(v)) c += (s - t) + v; else c += (v - t) + s;
    s = t;
  }
  return s + c;
}

template <class T>
T sum_squares(std::span<const T> x) { return dot_neumaier(x, x); }

template <class T>
T mean(std::span<const T> x) { RAID_REQUIRE(!x.empty(), insufficient_samples, "mean of empty span"); return sum_neumaier(x) / static_cast<T>(x.size()); }

/// Welford one-pass mean/variance (ddof = 1 by default).
struct Welford {
  usize n{0}; f64 mean{0}; f64 m2{0};
  void add(f64 x) { ++n; const f64 d = x - mean; mean += d / static_cast<f64>(n); m2 += d * (x - mean); }
  f64 variance(int ddof = 1) const { return n > static_cast<usize>(ddof) ? m2 / static_cast<f64>(n - static_cast<usize>(ddof)) : std::numeric_limits<f64>::quiet_NaN(); }
  /// Chan–Golub–LeVeque parallel combine.
  void merge(const Welford& o) {
    if (o.n == 0) return; if (n == 0) { *this = o; return; }
    const f64 d = o.mean - mean; const usize nn = n + o.n;
    mean += d * static_cast<f64>(o.n) / static_cast<f64>(nn);
    m2 += o.m2 + d * d * static_cast<f64>(n) * static_cast<f64>(o.n) / static_cast<f64>(nn);
    n = nn;
  }
};

template <class T>
T variance(std::span<const T> x, int ddof = 1) { Welford w; for (T v : x) w.add(static_cast<f64>(v)); return static_cast<T>(w.variance(ddof)); }

struct Moments { usize n; f64 mean, variance, skewness, kurtosis /* excess */, min, max; };
/// Two-pass moments; skewness/kurtosis use the biased (population) normalisation g1, g2 − 3 (as scipy with bias=True).
Moments moments(std::span<const f64> x);

/// Autocovariance at lag k with the 1/n normalisation (as statsmodels acovf(adjusted=False)).
f64 autocovariance(std::span<const f64> x, usize k, bool demean = true);
std::vector<f64> autocorrelation(std::span<const f64> x, usize max_lag);

/// Maximum drawdown of a cumulative-return path (fractional, ≥ 0), and the indices of the peak and trough.
struct Drawdown { f64 max_drawdown; usize peak_index; usize trough_index; };
Drawdown max_drawdown(std::span<const f64> cumulative);

}  // namespace raid::stats