1// CPU reference backend: FP64, Neumaier-compensated per-series reductions, parallel over series.2#include "raid/backend/backend.hpp"3#include "raid/core/thread_pool.hpp"4#include "raid/core/error.hpp"5#include "raid/statistics/reductions.hpp"6#include "raid/statistics/rolling.hpp"7#include "raid/models/garch.hpp"8#include <chrono>9#include <thread>10#include <cmath>11#ifdef __APPLE__12#include <sys/sysctl.h>13#endif1415namespace raid::backend {1617namespace {18struct 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(); } };19void 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; }20}2122class CpuBackend final : public Backend {23 public:24 explicit CpuBackend(unsigned threads) : pool_(threads == 0 ? nullptr : std::make_unique<ThreadPool>(threads)) {}25 std::string name() const override { return "cpu"; }26 raid::Backend kind() const override { return raid::Backend::cpu; }27 bool available() const override { return true; }28 Device device() const override {29 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());30#ifdef __APPLE__31 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";32 i64 mem = 0; n = sizeof mem; if (sysctlbyname("hw.memsize", &mem, &n, nullptr, 0) == 0) d.memory_bytes = static_cast<u64>(mem);33#else34 d.name = "cpu";35#endif36 d.details = "threads=" + std::to_string(pool().size());37 return d;38 }39 ThreadPool& pool() const { return pool_ ? *pool_ : ThreadPool::global(); }4041 template <class F> std::vector<f64> per_series(const Batch& b, F f) {42 std::vector<f64> out(b.n_series());43 pool().parallel_for(b.n_series(), [&](usize i0, usize i1) { for (usize i = i0; i < i1; ++i) out[i] = f(b.series(i)); });44 return out;45 }4647 std::vector<f64> sum_squares(const Batch& b, KernelStats* st) override {48 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;49 }50 std::vector<f64> sum_abs_adjacent(const Batch& b, KernelStats* st) override {51 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; });52 fill(st, "sum_abs_adjacent", t, b.values.size_bytes(), r.size() * 8); return r;53 }54 void semivariances(const Batch& b, std::vector<f64>& neg, std::vector<f64>& pos, KernelStats* st) override {55 Timer t; neg.assign(b.n_series(), 0); pos.assign(b.n_series(), 0);56 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; } });57 fill(st, "semivariances", t, b.values.size_bytes(), neg.size() * 16);58 }59 std::vector<f64> sum_fourth(const Batch& b, KernelStats* st) override {60 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; });61 fill(st, "sum_fourth", t, b.values.size_bytes(), r.size() * 8); return r;62 }63 std::vector<f64> rolling_sum_squares(const Batch& b, usize window, KernelStats* st) override {64 Timer t; std::vector<f64> out(b.values.size());65 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]); } });66 fill(st, "rolling_sum_squares", t, b.values.size_bytes(), out.size() * 8); return out;67 }68 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 {69 Timer t; const usize ns = offsets.size() - 1; xtx.assign(ns * k * k, 0); xty.assign(ns * k, 0); yty.assign(ns, 0);70 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]);71 std::vector<f64> acc(k * k, 0), accc(k * k, 0), ay(k, 0), ayc(k, 0); f64 syy = 0, cyy = 0;72 for (usize r = 0; r < n; ++r) { const f64* xr = x.data() + (b0 + r) * k; const f64 yr = y[b0 + r];73 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; }74 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; }75 const f64 q = yr * yr, tt = syy + q; if (std::abs(syy) >= q) cyy += (syy - tt) + q; else cyy += (q - tt) + syy; syy = tt; }76 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; }77 for (usize a = 0; a < k; ++a) xty[s * k + a] = ay[a] + ayc[a];78 yty[s] = syy + cyy; } });79 fill(st, "regression_stats", t, x.size_bytes() + y.size_bytes(), (xtx.size() + xty.size() + yty.size()) * 8);80 }81 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 {82 Timer t; const usize ns = returns.n_series(); std::vector<f64> out(ns * n_params);83 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); } });84 fill(st, "garch_nll_batch", t, returns.values.size_bytes() + params.size_bytes(), out.size() * 8); return out;85 }86 private:87 std::unique_ptr<ThreadPool> pool_;88};8990std::shared_ptr<Backend> make_cpu_backend(unsigned threads) { return std::make_shared<CpuBackend>(threads); }9192std::shared_ptr<Backend> resolve(const ExecutionPlan& plan, std::string* reason) {93 auto why = [&](const std::string& s) { if (reason) *reason = s; };94 if (plan.backend == raid::Backend::cpu) { why("cpu requested"); return make_cpu_backend(plan.threads); }95#ifdef RAID_ENABLE_METAL96 auto m = make_metal_backend(plan.precision);97 if (m && m->available()) { why(plan.backend == raid::Backend::metal ? "metal requested" : "auto: metal available"); return m; }98 if (plan.backend == raid::Backend::metal && !plan.allow_fallback) throw RaidError(ErrorCode::unsupported, "Metal backend requested but unavailable and fallback disabled");99 why("metal unavailable: falling back to cpu");100#else101 if (plan.backend == raid::Backend::metal && !plan.allow_fallback) throw RaidError(ErrorCode::unsupported, "RAID built without Metal support");102 why(plan.backend == raid::Backend::metal ? "built without Metal: falling back to cpu" : "auto: cpu (no Metal in this build)");103#endif104 return make_cpu_backend(plan.threads);105}106107std::vector<Device> list_devices() {108 std::vector<Device> d; d.push_back(make_cpu_backend()->device());109#ifdef RAID_ENABLE_METAL110 if (auto m = make_metal_backend()) if (m->available()) d.push_back(m->device());111#endif112 return d;113}114115} // namespace raid::backend116