// CPU reference backend: FP64, Neumaier-compensated per-series reductions, parallel over series.
#include "raid/backend/backend.hpp"
#include "raid/core/thread_pool.hpp"
#include "raid/core/error.hpp"
#include "raid/statistics/reductions.hpp"
#include "raid/statistics/rolling.hpp"
#include "raid/models/garch.hpp"
#include <chrono>
#include <thread>
#include <cmath>
#ifdef __APPLE__
#include <sys/sysctl.h>
#endif

namespace raid::backend {

namespace {
struct Timer { std::chrono::steady_clock::time_point t0 = std::chrono::steady_clock::now(); f64 s() const { return std::chrono::duration<f64>(std::chrono::steady_clock::now() - t0).count(); } };
void fill(KernelStats* st, const char* k, const Timer& t, usize bin, usize bout) { if (!st) return; st->kernel = k; st->backend = raid::Backend::cpu; st->precision = Precision::fp64; st->wall_seconds = t.s(); st->bytes_in = bin; st->bytes_out = bout; }
}

class CpuBackend final : public Backend {
 public:
  explicit CpuBackend(unsigned threads) : pool_(threads == 0 ? nullptr : std::make_unique<ThreadPool>(threads)) {}
  std::string name() const override { return "cpu"; }
  raid::Backend kind() const override { return raid::Backend::cpu; }
  bool available() const override { return true; }
  Device device() const override {
    Device d; d.kind = raid::Backend::cpu; d.supports_fp64 = true; d.unified_memory = true; d.compute_units = static_cast<int>(std::thread::hardware_concurrency());
#ifdef __APPLE__
    char buf[256]; size_t n = sizeof buf; if (sysctlbyname("machdep.cpu.brand_string", buf, &n, nullptr, 0) == 0) d.name = std::string(buf, n - 1); else d.name = "cpu";
    i64 mem = 0; n = sizeof mem; if (sysctlbyname("hw.memsize", &mem, &n, nullptr, 0) == 0) d.memory_bytes = static_cast<u64>(mem);
#else
    d.name = "cpu";
#endif
    d.details = "threads=" + std::to_string(pool().size());
    return d;
  }
  ThreadPool& pool() const { return pool_ ? *pool_ : ThreadPool::global(); }

  template <class F> std::vector<f64> per_series(const Batch& b, F f) {
    std::vector<f64> out(b.n_series());
    pool().parallel_for(b.n_series(), [&](usize i0, usize i1) { for (usize i = i0; i < i1; ++i) out[i] = f(b.series(i)); });
    return out;
  }

  std::vector<f64> sum_squares(const Batch& b, KernelStats* st) override {
    Timer t; auto r = per_series(b, [](std::span<const f64> x) { return stats::sum_squares(x); }); fill(st, "sum_squares", t, b.values.size_bytes(), r.size() * 8); return r;
  }
  std::vector<f64> sum_abs_adjacent(const Batch& b, KernelStats* st) override {
    Timer t; auto r = per_series(b, [](std::span<const f64> x) { f64 s = 0, c = 0; for (usize i = 1; i < x.size(); ++i) { const f64 v = std::abs(x[i]) * std::abs(x[i - 1]), tt = s + v; if (std::abs(s) >= std::abs(v)) c += (s - tt) + v; else c += (v - tt) + s; s = tt; } return s + c; });
    fill(st, "sum_abs_adjacent", t, b.values.size_bytes(), r.size() * 8); return r;
  }
  void semivariances(const Batch& b, std::vector<f64>& neg, std::vector<f64>& pos, KernelStats* st) override {
    Timer t; neg.assign(b.n_series(), 0); pos.assign(b.n_series(), 0);
    pool().parallel_for(b.n_series(), [&](usize i0, usize i1) { for (usize i = i0; i < i1; ++i) { f64 sn = 0, cn = 0, sp = 0, cp = 0; for (f64 v : b.series(i)) { const f64 q = v * v; if (v < 0) { const f64 tt = sn + q; if (std::abs(sn) >= q) cn += (sn - tt) + q; else cn += (q - tt) + sn; sn = tt; } else if (v > 0) { const f64 tt = sp + q; if (std::abs(sp) >= q) cp += (sp - tt) + q; else cp += (q - tt) + sp; sp = tt; } } neg[i] = sn + cn; pos[i] = sp + cp; } });
    fill(st, "semivariances", t, b.values.size_bytes(), neg.size() * 16);
  }
  std::vector<f64> sum_fourth(const Batch& b, KernelStats* st) override {
    Timer t; auto r = per_series(b, [](std::span<const f64> x) { f64 s = 0, c = 0; for (f64 v : x) { const f64 q = v * v * v * v, tt = s + q; if (std::abs(s) >= q) c += (s - tt) + q; else c += (q - tt) + s; s = tt; } return s + c; });
    fill(st, "sum_fourth", t, b.values.size_bytes(), r.size() * 8); return r;
  }
  std::vector<f64> rolling_sum_squares(const Batch& b, usize window, KernelStats* st) override {
    Timer t; std::vector<f64> out(b.values.size());
    pool().parallel_for(b.n_series(), [&](usize i0, usize i1) { for (usize i = i0; i < i1; ++i) { auto x = b.series(i); std::vector<f64> sq(x.size()); for (usize j = 0; j < x.size(); ++j) sq[j] = x[j] * x[j]; auto r = stats::rolling_sum(sq, stats::RollingOptions{window, 0, 1024}); std::copy(r.begin(), r.end(), out.begin() + b.offsets[i]); } });
    fill(st, "rolling_sum_squares", t, b.values.size_bytes(), out.size() * 8); return out;
  }
  void regression_stats(std::span<const f64> x, std::span<const f64> y, std::span<const i64> offsets, usize k, std::vector<f64>& xtx, std::vector<f64>& xty, std::vector<f64>& yty, KernelStats* st) override {
    Timer t; const usize ns = offsets.size() - 1; xtx.assign(ns * k * k, 0); xty.assign(ns * k, 0); yty.assign(ns, 0);
    pool().parallel_for(ns, [&](usize i0, usize i1) { for (usize s = i0; s < i1; ++s) { const usize b0 = static_cast<usize>(offsets[s]), n = static_cast<usize>(offsets[s + 1] - offsets[s]);
        std::vector<f64> acc(k * k, 0), accc(k * k, 0), ay(k, 0), ayc(k, 0); f64 syy = 0, cyy = 0;
        for (usize r = 0; r < n; ++r) { const f64* xr = x.data() + (b0 + r) * k; const f64 yr = y[b0 + r];
          for (usize a = 0; a < k; ++a) { for (usize c = a; c < k; ++c) { const f64 v = xr[a] * xr[c], tt = acc[a * k + c] + v; if (std::abs(acc[a * k + c]) >= std::abs(v)) accc[a * k + c] += (acc[a * k + c] - tt) + v; else accc[a * k + c] += (v - tt) + acc[a * k + c]; acc[a * k + c] = tt; }
            const f64 v = xr[a] * yr, tt = ay[a] + v; if (std::abs(ay[a]) >= std::abs(v)) ayc[a] += (ay[a] - tt) + v; else ayc[a] += (v - tt) + ay[a]; ay[a] = tt; }
          const f64 q = yr * yr, tt = syy + q; if (std::abs(syy) >= q) cyy += (syy - tt) + q; else cyy += (q - tt) + syy; syy = tt; }
        for (usize a = 0; a < k; ++a) for (usize c = a; c < k; ++c) { const f64 v = acc[a * k + c] + accc[a * k + c]; xtx[s * k * k + a * k + c] = v; xtx[s * k * k + c * k + a] = v; }
        for (usize a = 0; a < k; ++a) xty[s * k + a] = ay[a] + ayc[a];
        yty[s] = syy + cyy; } });
    fill(st, "regression_stats", t, x.size_bytes() + y.size_bytes(), (xtx.size() + xty.size() + yty.size()) * 8);
  }
  std::vector<f64> garch_nll_batch(const Batch& returns, std::span<const f64> params, usize n_params, usize p, int model_code, f64 nu_fixed, KernelStats* st) override {
    Timer t; const usize ns = returns.n_series(); std::vector<f64> out(ns * n_params);
    pool().parallel_for(ns * n_params, [&](usize i0, usize i1) { for (usize idx = i0; idx < i1; ++idx) { const usize s = idx / n_params, q = idx % n_params; const f64* th = params.data() + (s * n_params + q) * p; out[idx] = models::garch_nll(returns.series(s), std::span<const f64>(th, p), static_cast<models::GarchModel>(model_code), nu_fixed); } });
    fill(st, "garch_nll_batch", t, returns.values.size_bytes() + params.size_bytes(), out.size() * 8); return out;
  }
 private:
  std::unique_ptr<ThreadPool> pool_;
};

std::shared_ptr<Backend> make_cpu_backend(unsigned threads) { return std::make_shared<CpuBackend>(threads); }

std::shared_ptr<Backend> resolve(const ExecutionPlan& plan, std::string* reason) {
  auto why = [&](const std::string& s) { if (reason) *reason = s; };
  if (plan.backend == raid::Backend::cpu) { why("cpu requested"); return make_cpu_backend(plan.threads); }
#ifdef RAID_ENABLE_METAL
  auto m = make_metal_backend(plan.precision);
  if (m && m->available()) { why(plan.backend == raid::Backend::metal ? "metal requested" : "auto: metal available"); return m; }
  if (plan.backend == raid::Backend::metal && !plan.allow_fallback) throw RaidError(ErrorCode::unsupported, "Metal backend requested but unavailable and fallback disabled");
  why("metal unavailable: falling back to cpu");
#else
  if (plan.backend == raid::Backend::metal && !plan.allow_fallback) throw RaidError(ErrorCode::unsupported, "RAID built without Metal support");
  why(plan.backend == raid::Backend::metal ? "built without Metal: falling back to cpu" : "auto: cpu (no Metal in this build)");
#endif
  return make_cpu_backend(plan.threads);
}

std::vector<Device> list_devices() {
  std::vector<Device> d; d.push_back(make_cpu_backend()->device());
#ifdef RAID_ENABLE_METAL
  if (auto m = make_metal_backend()) if (m->available()) d.push_back(m->device());
#endif
  return d;
}

}  // namespace raid::backend