#include "raid/statistics/reductions.hpp"

namespace raid::stats {

Moments moments(std::span<const f64> x) {
  RAID_REQUIRE(x.size() >= 2, insufficient_samples, "moments need n >= 2");
  const f64 n = static_cast<f64>(x.size());
  const f64 m = sum_neumaier(x) / n;
  std::vector<f64> d2(x.size()), d3(x.size()), d4(x.size());
  f64 mn = x[0], mx = x[0];
  for (usize i = 0; i < x.size(); ++i) {
    const f64 d = x[i] - m; d2[i] = d * d; d3[i] = d2[i] * d; d4[i] = d2[i] * d2[i];
    mn = std::min(mn, x[i]); mx = std::max(mx, x[i]);
  }
  const f64 m2 = sum_neumaier(std::span<const f64>(d2)) / n, m3 = sum_neumaier(std::span<const f64>(d3)) / n, m4 = sum_neumaier(std::span<const f64>(d4)) / n;
  const f64 var_unb = sum_neumaier(std::span<const f64>(d2)) / (n - 1.0);
  return {x.size(), m, var_unb, m2 > 0 ? m3 / std::pow(m2, 1.5) : std::numeric_limits<f64>::quiet_NaN(), m2 > 0 ? m4 / (m2 * m2) - 3.0 : std::numeric_limits<f64>::quiet_NaN(), mn, mx};
}

f64 autocovariance(std::span<const f64> x, usize k, bool demean) {
  RAID_REQUIRE(k < x.size(), insufficient_samples, "autocovariance lag >= n");
  const f64 n = static_cast<f64>(x.size());
  const f64 m = demean ? sum_neumaier(x) / n : 0.0;
  std::vector<f64> p(x.size() - k);
  for (usize i = 0; i + k < x.size(); ++i) p[i] = (x[i] - m) * (x[i + k] - m);
  return sum_neumaier(std::span<const f64>(p)) / n;
}

std::vector<f64> autocorrelation(std::span<const f64> x, usize max_lag) {
  const f64 c0 = autocovariance(x, 0);
  std::vector<f64> r(max_lag + 1, 1.0);
  for (usize k = 1; k <= max_lag && k < x.size(); ++k) r[k] = c0 > 0 ? autocovariance(x, k) / c0 : std::numeric_limits<f64>::quiet_NaN();
  return r;
}

Drawdown max_drawdown(std::span<const f64> cum) {
  Drawdown d{0.0, 0, 0};
  if (cum.empty()) return d;
  f64 peak = cum[0]; usize peak_i = 0;
  for (usize i = 0; i < cum.size(); ++i) {
    if (cum[i] > peak) { peak = cum[i]; peak_i = i; }
    const f64 dd = peak != 0.0 ? (peak - cum[i]) / peak : 0.0;
    if (dd > d.max_drawdown) { d.max_drawdown = dd; d.peak_index = peak_i; d.trough_index = i; }
  }
  return d;
}

}  // namespace raid::stats