#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <nlohmann/json.hpp>
#include <fstream>
#include "raid/estimators/har.hpp"
#include "raid/core/rng.hpp"
using namespace raid; using json = nlohmann::json; using Catch::Matchers::WithinRel; using Catch::Matchers::WithinAbs;

static void load_rv(std::vector<double>& rv, std::vector<double>& rq) { std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/har_rv.csv"); REQUIRE(f.good()); std::string line; std::getline(f, line); while (std::getline(f, line)) { auto c = line.find(','); rv.push_back(std::stod(line.substr(0, c))); rq.push_back(std::stod(line.substr(c + 1))); } }

TEST_CASE("HAR-RV / log-HAR / HARQ with HAC match statsmodels on the fixture", "[reference][har]") {
  std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/reference_values.json"); json ref; f >> ref; auto har = ref["har"];
  std::vector<double> rv, rq; load_rv(rv, rq);
  struct C { const char* key; int h; bool log; bool harq; }; const C cases[] = {{"h1", 1, false, false}, {"h5", 5, false, false}, {"h22", 22, false, false}, {"log_h1", 1, true, false}, {"harq_h1", 1, false, true}};
  for (auto& c : cases) {
    est::HarSpec s; s.horizon = c.h; s.log = c.log; s.harq = c.harq; auto e = est::har_fit(rv, s, c.harq ? std::span<const double>(rq) : std::span<const double>());
    INFO(c.key); REQUIRE(e.ok()); auto r = har[c.key]; REQUIRE(e.n_obs == r["n"].get<usize>());
    for (usize j = 0; j < e.coefficients.size(); ++j) { CHECK_THAT(e.coefficients[j], WithinRel(r["params"][j].get<double>(), 1e-9)); CHECK_THAT(e.inference->std_errors[j], WithinRel(r["hac_bse"][j].get<double>(), 1e-8)); }
    CHECK(e.inference->covariance_type == "hac_bartlett(" + std::to_string(r["hac_lags"].get<int>()) + ")"); CHECK_THAT(e.loss, WithinRel(r["ssr"].get<double>(), 1e-9));
  }
}
TEST_CASE("HAR batched (sufficient statistics, auto backend) equals per-series QR fits", "[har]") {
  std::vector<double> rv, rq; load_rv(rv, rq); const usize S = 6; std::vector<double> cat; std::vector<i64> off = {0};
  Rng r(5); for (usize s = 0; s < S; ++s) { for (double v : rv) cat.push_back(v * (1 + 0.1 * r.normal()) * (s + 1)); off.push_back(static_cast<i64>(cat.size())); }
  est::HarSpec spec; spec.horizon = 1; spec.cov.type = stats::CovType::classical;
  for (auto bk : {Backend::cpu, Backend::automatic}) { ExecutionPlan p; p.backend = bk; auto b = est::har_batch(cat, off, spec, p); REQUIRE(b.fits.size() == S);
    for (usize s = 0; s < S; ++s) { auto single = est::har_fit(std::span<const double>(cat).subspan(static_cast<usize>(off[s]), static_cast<usize>(off[s + 1] - off[s])), spec); REQUIRE(b.fits[s].ok());
      for (usize j = 0; j < 4; ++j) { INFO("backend " << to_string(bk) << " s=" << s << " j=" << j); CHECK_THAT(b.fits[s].coefficients[j], WithinAbs(single.coefficients[j], bk == Backend::cpu ? 1e-9 : 1e-4 * std::max(1.0, std::abs(single.coefficients[j])))); } } }
}
TEST_CASE("HAR rolling forecasts are causal and losses are finite", "[har]") {
  std::vector<double> rv, rq; load_rv(rv, rq); est::HarSpec s; s.horizon = 5; auto roll = est::har_rolling(rv, s, 500, 10);
  CHECK(roll.n_fits > 50); for (usize i = 0; i < roll.origin.size(); ++i) { CHECK(roll.forecast[i] > 0); double a = 0; for (int j = 1; j <= 5; ++j) a += rv[roll.origin[i] + j]; CHECK_THAT(roll.actual[i], WithinRel(a / 5.0, 1e-12)); }
  auto l = est::forecast_loss(roll.actual, roll.forecast); CHECK(l.n == roll.actual.size()); CHECK(std::isfinite(l.qlike)); CHECK(l.qlike >= 0);
  // the forecast at origin t uses only rv[0..t]: refitting with the future perturbed must not change it
  std::vector<double> rv2 = rv; for (usize t = roll.origin[0] + 1; t < rv2.size(); ++t) rv2[t] *= 3.0; auto roll2 = est::har_rolling(rv2, s, 500, 10); CHECK_THAT(roll2.forecast[0], WithinRel(roll.forecast[0], 1e-12));
  est::HarSpec lg; lg.log = true; auto e = est::har_fit(rv, lg); CHECK(e.extra.contains("smearing")); CHECK(est::har_forecast(e, rv, lg) > 0);
}