// E5 — inference at scale and reproducibility: (a) GARCH(1,1) Monte Carlo R = 1 000 (CPU parallel): bias, RMSE, coverage,
// failures, time; (b) HAR-RV h=1/5 HAC coverage R = 1 000; (c) stationary-bootstrap RV intervals on real sessions
// (coverage of the "long-sample" IV proxy is not identifiable — we report interval widths and reproducibility);
// (d) bitwise reproducibility: same seed → identical bootstrap statistics across thread counts, and identical Metal
// realized measures across runs; CPU-vs-Metal differences reported in ulps-equivalent relative error.
// Usage: bench_e5_inference <out.jsonl> [R]
#include <cstdio>
#include <fstream>
#include <atomic>
#include <nlohmann/json.hpp>
#include "raid/core/manifest.hpp"
#include "raid/core/rng.hpp"
#include "raid/core/thread_pool.hpp"
#include "raid/models/garch.hpp"
#include "raid/estimators/har.hpp"
#include "raid/estimators/realized.hpp"
#include "raid/statistics/bootstrap.hpp"
#include "raid/statistics/reductions.hpp"
using namespace raid; using json = nlohmann::json;
namespace { f64 now() { return std::chrono::duration<f64>(std::chrono::steady_clock::now().time_since_epoch()).count(); } }

int main(int argc, char** argv) {
  if (argc < 2) { std::fprintf(stderr, "usage: %s out.jsonl [R]\n", argv[0]); return 2; }
  const int R = argc > 2 ? std::atoi(argv[2]) : 1000; std::ofstream out(argv[1], std::ios::app);
  ExperimentManifest man; man.experiment = "E5.inference"; man.created_at_utc = now_utc_iso(); man.raid_version = RAID_VERSION; man.host = host_info(); man.seed = 20261102; man.parameters = {{"R", std::to_string(R)}};
  // (a) GARCH MC
  for (auto& [n, label] : std::vector<std::pair<usize, const char*>>{{500, "n500"}, {2000, "n2000"}}) { const double tr[4] = {0.0, 0.02, 0.08, 0.90}; std::vector<std::array<double, 4>> est(R), se(R); std::vector<int> okv(R, 0), its(R, 0); const f64 t0 = now();
    ThreadPool::global().parallel_for(R, [&](usize a, usize b) { for (usize rep = a; rep < b; ++rep) { Rng rng(20261102, rep); std::vector<double> r(n); double s2 = tr[1] / (1 - tr[2] - tr[3]), e = 0; for (usize t = 0; t < n; ++t) { s2 = tr[1] + tr[2] * e * e + tr[3] * s2; e = std::sqrt(s2) * rng.normal(); r[t] = e; } models::GarchFitOptions o; auto f = models::garch_fit(r, o); its[rep] = f.convergence.iterations; if (!f.ok()) continue; okv[rep] = 1; for (int j = 0; j < 4; ++j) { est[rep][j] = f.coefficients[j]; se[rep][j] = f.inference->std_errors[j]; } } });
    json row = {{"kind", "garch_mc"}, {"experiment", "E5"}, {"n", n}, {"R", R}, {"wall_seconds", now() - t0}, {"threads", ThreadPool::global().size()}}; int ok = 0; for (int v : okv) ok += v; row["converged"] = ok;
    for (int j = 0; j < 4; ++j) { double bias = 0, mse = 0; int cov = 0; for (int rep = 0; rep < R; ++rep) { if (!okv[rep]) continue; double d = est[rep][j] - tr[j]; bias += d; mse += d * d; if (std::abs(d) <= 1.96 * se[rep][j]) ++cov; } const char* nm[4] = {"mu", "omega", "alpha", "beta"}; row[nm[j]] = {{"bias", bias / ok}, {"rmse", std::sqrt(mse / ok)}, {"coverage95", static_cast<double>(cov) / ok}}; }
    out << row.dump() << "\n"; out.flush(); std::fprintf(stderr, "E5 garch MC %s: %d/%d converged %.1fs\n", label, ok, R, now() - t0); }
  // (b) HAR HAC coverage
  for (int h : {1, 5}) { auto simulate = [](usize nn, u64 seed) { Rng g(4242, seed); std::vector<double> rv(nn); double lv = std::log(1e-4); for (usize t = 0; t < nn; ++t) { lv = 0.3 * std::log(1e-4) + 0.7 * lv + 0.4 * g.normal(); rv[t] = std::exp(lv); } return rv; };
    est::HarSpec spec; spec.horizon = h; auto big = simulate(400000, 999); auto pseudo = est::har_fit(big, spec); std::vector<std::array<int, 4>> cov(R); const f64 t0 = now();
    ThreadPool::global().parallel_for(R, [&](usize a, usize b) { for (usize rep = a; rep < b; ++rep) { auto rv = simulate(1500, rep); auto f = est::har_fit(rv, spec); for (int j = 0; j < 4; ++j) cov[rep][j] = std::abs(f.coefficients[j] - pseudo.coefficients[j]) <= 1.96 * f.inference->std_errors[j]; } });
    json row = {{"kind", "har_mc"}, {"experiment", "E5"}, {"h", h}, {"T", 1500}, {"R", R}, {"wall_seconds", now() - t0}}; for (int j = 0; j < 4; ++j) { int c = 0; for (int rep = 0; rep < R; ++rep) c += cov[rep][j]; row["coverage95"][j] = static_cast<double>(c) / R; } out << row.dump() << "\n"; out.flush(); std::fprintf(stderr, "E5 har h=%d done\n", h); }
  // (d) reproducibility
  { Rng g(99); std::vector<double> x(5000); for (auto& v : x) v = 1e-3 * g.student_t(5.0); auto st = [](std::span<const double> v) { return stats::sum_squares<double>(v); };
    stats::BootstrapOptions o; o.replicas = 2000; o.mean_block = 10; std::vector<std::vector<double>> results; std::vector<unsigned> thr = {1, 4, 18}; for (unsigned t : thr) { o.threads = t; auto b = stats::bootstrap(x, st, o); results.push_back(b.statistics); }
    bool same = results[0] == results[1] && results[1] == results[2]; out << json({{"kind", "repro_bootstrap"}, {"experiment", "E5"}, {"replicas", 2000}, {"threads", thr}, {"bitwise_identical_across_threads", same}}).dump() << "\n";
    std::vector<i64> off; for (usize s = 0; s <= 50; ++s) off.push_back(static_cast<i64>(s * 100)); backend::Batch b{x, off}; est::RealizedOptions rc; rc.plan.backend = Backend::cpu; auto c1 = est::realized_measures_batch(b, rc); rc.plan.threads = 1; auto c2 = est::realized_measures_batch(b, rc);
    json rr = {{"kind", "repro_realized"}, {"experiment", "E5"}, {"cpu_threads_vs_1_identical", c1.rv == c2.rv}};
#ifdef RAID_ENABLE_METAL
    est::RealizedOptions rm; rm.plan.backend = Backend::metal; auto m1 = est::realized_measures_batch(b, rm), m2 = est::realized_measures_batch(b, rm); rr["metal_runs_identical"] = m1.rv == m2.rv && m1.bpv == m2.bpv && m1.rq == m2.rq; double maxrel = 0; for (usize s = 0; s < c1.rv.size(); ++s) maxrel = std::max(maxrel, std::abs(m1.rv[s] - c1.rv[s]) / c1.rv[s]); rr["metal_vs_cpu_max_rel_rv"] = maxrel; rr["metal_backend_used"] = std::string(to_string(m1.stats.backend));
#endif
    out << rr.dump() << "\n"; }
  man.save(std::string(argv[1]) + ".manifest.json"); return 0;
}