SPB Git
spboucher.ai

spb/wp40_uqo

Public MIT

WP40 RAID: Phase 3 — Metal backend (double-float accumulation, fused realized kernel, rolling sums, regression stats, batched GARCH NLL), accuracy/determinism tests, primitives microbenchmark + analysis script

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Simon-Pierre Boucher committed (Nov 2, 2026) parent 516afa5 Browse files

9 changed files +607 −8

9 changed files
Unified Split
modified CMakeLists.txt +1 −1
@@ -112,7 +112,7 @@ if(RAID_ENABLE_METAL)
112 foreach(s ${RAID_METAL_SHADERS}) 112 foreach(s ${RAID_METAL_SHADERS})
113 get_filename_component(n ${s} NAME_WE) 113 get_filename_component(n ${s} NAME_WE)
114 add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${n}.air 114 add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${n}.air
115− COMMAND ${RAID_XCRUN} -sdk macosx metal -std=metal3.2 -O2 -c ${s} -o ${CMAKE_CURRENT_BINARY_DIR}/${n}.air 115+ COMMAND ${RAID_XCRUN} -sdk macosx metal -O2 -fno-fast-math -c ${s} -o ${CMAKE_CURRENT_BINARY_DIR}/${n}.air
116 DEPENDS ${s} COMMENT "metal ${n}") 116 DEPENDS ${s} COMMENT "metal ${n}")
117 list(APPEND _airs ${CMAKE_CURRENT_BINARY_DIR}/${n}.air) 117 list(APPEND _airs ${CMAKE_CURRENT_BINARY_DIR}/${n}.air)
118 endforeach() 118 endforeach()
added backends/metal/host/metal_backend.mm +241 −0
@@ -0,0 +1,241 @@
1+// Metal backend host code (Objective-C++). See shaders/raid_kernels.metal for the kernels and the precision design.
2+#import <Metal/Metal.h>
3+#import <Foundation/Foundation.h>
4+#include "raid/backend/backend.hpp"
5+#include "raid/core/error.hpp"
6+#include "raid/core/log.hpp"
7+#include "raid/statistics/reductions.hpp"
8+#include <chrono>
9+#include <cmath>
10+#include <map>
11+#include <mutex>
12+#include <unistd.h>
13+#include <mach-o/dyld.h>
14+
15+namespace raid { namespace metal { const char* msl_bundle(); } }
16+
17+namespace raid::backend {
18+
19+namespace {
20+struct Chunk { uint32_t seg, begin, end, pad; };
21+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(); } };
22+}
23+
24+class MetalBackend final : public Backend {
25+ public:
26+ explicit MetalBackend(PrecisionPolicy policy) : policy_(policy) {
27+ @autoreleasepool {
28+ device_ = MTLCreateSystemDefaultDevice();
29+ if (!device_) { log::warn("metal: no device"); return; }
30+ queue_ = [device_ newCommandQueue];
31+ available_ = device_ != nil && queue_ != nil;
32+ }
33+ }
34+ ~MetalBackend() override = default;
35+ std::string name() const override { return "metal"; }
36+ raid::Backend kind() const override { return raid::Backend::metal; }
37+ bool available() const override { return available_; }
38+ Device device() const override {
39+ Device d; d.kind = raid::Backend::metal; d.supports_fp64 = false; d.unified_memory = device_.hasUnifiedMemory;
40+ d.name = device_ ? std::string([device_.name UTF8String]) : "none"; d.memory_bytes = device_.recommendedMaxWorkingSetSize; d.max_buffer_bytes = device_.maxBufferLength;
41+ d.compute_units = 0; // Metal does not expose the GPU core count; recorded from system_profiler in the manifest
42+ d.details = std::string("families:") + ([device_ supportsFamily:MTLGPUFamilyApple9] ? "Apple9" : ([device_ supportsFamily:MTLGPUFamilyApple8] ? "Apple8" : "Apple7-")) + " metal3:" + ([device_ supportsFamily:MTLGPUFamilyMetal3] ? "yes" : "no") + " library:" + library_source_;
43+ return d;
44+ }
45+
46+ // ---- library / pipelines ----------------------------------------------------------------------------------------
47+ bool ensure_library(KernelStats* st) {
48+ std::lock_guard<std::mutex> lk(m_);
49+ if (library_) return true;
50+ @autoreleasepool {
51+ NSError* err = nil; Timer t;
52+ // 1. prebuilt metallib: RAID_METALLIB env, compile-time path, next to the executable
53+ std::vector<std::string> candidates;
54+ if (const char* e = std::getenv("RAID_METALLIB")) candidates.push_back(e);
55+#ifdef RAID_METALLIB_PATH
56+ candidates.push_back(RAID_METALLIB_PATH);
57+#endif
58+ { char buf[4096]; uint32_t sz = sizeof buf; if (_NSGetExecutablePath(buf, &sz) == 0) { std::string p(buf); auto k = p.rfind('/'); if (k != std::string::npos) candidates.push_back(p.substr(0, k) + "/raid.metallib"); } }
59+ for (auto& c : candidates) {
60+ if (access(c.c_str(), R_OK) != 0) continue;
61+ library_ = [device_ newLibraryWithURL:[NSURL fileURLWithPath:[NSString stringWithUTF8String:c.c_str()]] error:&err];
62+ if (library_) { library_source_ = "metallib:" + c; compile_seconds_ = t.s(); break; }
63+ }
64+ // 2. runtime compilation of the embedded MSL
65+ if (!library_) {
66+ MTLCompileOptions* opts = [MTLCompileOptions new]; opts.fastMathEnabled = NO; // keep IEEE semantics (no reassociation): required by the double-float transforms
67+ library_ = [device_ newLibraryWithSource:[NSString stringWithUTF8String:raid::metal::msl_bundle()] options:opts error:&err];
68+ if (!library_) { log::error("metal: shader compilation failed", {{"error", err ? std::string([[err localizedDescription] UTF8String]) : "?"}}); available_ = false; return false; }
69+ library_source_ = "runtime"; compile_seconds_ = t.s(); if (st) st->compiled_now = true;
70+ }
71+ }
72+ return true;
73+ }
74+ id<MTLComputePipelineState> pipeline(const std::string& fn) {
75+ std::lock_guard<std::mutex> lk(m_);
76+ auto it = pipelines_.find(fn); if (it != pipelines_.end()) return it->second;
77+ @autoreleasepool {
78+ NSError* err = nil; id<MTLFunction> f = [library_ newFunctionWithName:[NSString stringWithUTF8String:fn.c_str()]];
79+ RAID_REQUIRE(f != nil, unsupported, "metal: kernel not found: " + fn);
80+ id<MTLComputePipelineState> p = [device_ newComputePipelineStateWithFunction:f error:&err];
81+ RAID_REQUIRE(p != nil, unsupported, "metal: pipeline failed for " + fn + ": " + (err ? std::string([[err localizedDescription] UTF8String]) : "?"));
82+ pipelines_[fn] = p; return p;
83+ }
84+ }
85+ // ---- buffer pool (shared storage; unified memory) ---------------------------------------------------------------
86+ id<MTLBuffer> buffer(usize bytes, usize* allocs) {
87+ std::lock_guard<std::mutex> lk(m_);
88+ bytes = std::max<usize>(bytes, 16); usize cls = 1; while (cls < bytes) cls <<= 1;
89+ auto& v = pool_[cls]; if (!v.empty()) { auto b = v.back(); v.pop_back(); return b; }
90+ if (allocs) ++*allocs; ++total_allocs_;
91+ id<MTLBuffer> b = [device_ newBufferWithLength:cls options:MTLResourceStorageModeShared];
92+ RAID_REQUIRE(b != nil, memory, "metal: buffer allocation failed (" + std::to_string(cls) + " bytes)");
93+ return b;
94+ }
95+ void release(id<MTLBuffer> b) { std::lock_guard<std::mutex> lk(m_); pool_[b.length].push_back(b); }
96+ usize total_allocations() const { return total_allocs_; }
97+
98+ // ---- helpers -----------------------------------------------------------------------------------------------------
99+ std::vector<Chunk> make_chunks(std::span<const i64> offsets, usize chunk_max, usize* n_segments) const {
100+ std::vector<Chunk> ch; const usize S = offsets.size() - 1; if (n_segments) *n_segments = S;
101+ for (usize s = 0; s < S; ++s) { const u64 b = static_cast<u64>(offsets[s]), e = static_cast<u64>(offsets[s + 1]); RAID_REQUIRE(e <= 0xFFFFFFFFull, unsupported, "metal: segments beyond 2^32 elements"); if (e == b) { ch.push_back({static_cast<uint32_t>(s), static_cast<uint32_t>(b), static_cast<uint32_t>(e), 0}); continue; }
102+ for (u64 x = b; x < e; x += chunk_max) ch.push_back({static_cast<uint32_t>(s), static_cast<uint32_t>(x), static_cast<uint32_t>(std::min<u64>(e, x + chunk_max)), 0}); }
103+ return ch;
104+ }
105+ usize chunk_size() const { return policy_.fp32_chunk_max == 0 ? 4096 : policy_.fp32_chunk_max; }
106+ static void finish(KernelStats* st, const char* k, const Timer& t, id<MTLCommandBuffer> cb, usize bin, usize bout, usize allocs, const std::string& note) {
107+ if (!st) return; st->kernel = k; st->backend = raid::Backend::metal; st->precision = Precision::fp32; st->wall_seconds = t.s(); st->device_seconds = cb ? (cb.GPUEndTime - cb.GPUStartTime) : 0; st->bytes_in = bin; st->bytes_out = bout; st->buffer_allocations = allocs; st->note = note;
108+ }
109+ // Fused realized kernel; returns 5 vectors combined in FP64 on the CPU (adds the chunk-boundary |x_i||x_{i-1}| terms).
110+ void realized_all(const Batch& b, std::vector<f64>& s2, std::vector<f64>& adj, std::vector<f64>& neg, std::vector<f64>& pos, std::vector<f64>& s4, KernelStats* st, int single_op = -1) {
111+ Timer t; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
112+ usize S; auto chunks = make_chunks(b.offsets, chunk_size(), &S); const usize n = b.values.size(), C = chunks.size();
113+ s2.assign(S, 0); adj.assign(S, 0); neg.assign(S, 0); pos.assign(S, 0); s4.assign(S, 0);
114+ if (n == 0 || C == 0) { finish(st, "realized_fused", t, nil, 0, 0, 0, library_source_); return; }
115+ @autoreleasepool {
116+ id<MTLBuffer> bx = buffer(n * 4, &allocs), bc = buffer(C * sizeof(Chunk), &allocs), bo = buffer(C * 5 * 8, &allocs);
117+ Timer tconv; float* xf = static_cast<float*>(bx.contents); for (usize i = 0; i < n; ++i) xf[i] = static_cast<float>(b.values[i]); const f64 conv_s = tconv.s();
118+ std::memcpy(bc.contents, chunks.data(), C * sizeof(Chunk));
119+ id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder];
120+ const uint32_t op = single_op < 0 ? 0 : static_cast<uint32_t>(single_op);
121+ id<MTLComputePipelineState> p = pipeline(single_op < 0 ? "realized_fused" : "segment_reduce");
122+ [enc setComputePipelineState:p]; [enc setBuffer:bx offset:0 atIndex:0]; [enc setBuffer:bc offset:0 atIndex:1]; [enc setBuffer:bo offset:0 atIndex:2]; if (single_op >= 0) [enc setBytes:&op length:4 atIndex:3];
123+ [enc dispatchThreadgroups:MTLSizeMake(C, 1, 1) threadsPerThreadgroup:MTLSizeMake(256, 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
124+ RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
125+ const float* o = static_cast<const float*>(bo.contents);
126+ for (usize c = 0; c < C; ++c) { const usize s = chunks[c].seg;
127+ if (single_op < 0) { s2[s] += static_cast<f64>(o[c * 10 + 0]) + o[c * 10 + 1]; adj[s] += static_cast<f64>(o[c * 10 + 2]) + o[c * 10 + 3]; neg[s] += static_cast<f64>(o[c * 10 + 4]) + o[c * 10 + 5]; pos[s] += static_cast<f64>(o[c * 10 + 6]) + o[c * 10 + 7]; s4[s] += static_cast<f64>(o[c * 10 + 8]) + o[c * 10 + 9]; }
128+ else { f64 v = static_cast<f64>(o[c * 2]) + o[c * 2 + 1]; (single_op == 0 ? s2 : single_op == 1 ? adj : single_op == 2 ? neg : single_op == 3 ? pos : s4)[s] += v; }
129+ // chunk boundary term for the adjacent product (FP64 on the CPU)
130+ if ((single_op < 0 || single_op == 1) && chunks[c].begin > static_cast<uint32_t>(b.offsets[s])) adj[s] += std::abs(b.values[chunks[c].begin]) * std::abs(b.values[chunks[c].begin - 1]); }
131+ finish(st, single_op < 0 ? "realized_fused" : "segment_reduce", t, cb, n * 8, S * 40, allocs, library_source_ + " chunks=" + std::to_string(C) + " convert_s=" + std::to_string(conv_s));
132+ release(bx); release(bc); release(bo);
133+ }
134+ }
135+
136+ // ---- Backend interface ---------------------------------------------------------------------------------------------
137+ std::vector<f64> sum_squares(const Batch& b, KernelStats* st) override { std::vector<f64> a, c, d, e, f; realized_all(b, a, c, d, e, f, st, 0); return a; }
138+ std::vector<f64> sum_abs_adjacent(const Batch& b, KernelStats* st) override { std::vector<f64> a, c, d, e, f; realized_all(b, a, c, d, e, f, st, 1); return c; }
139+ void semivariances(const Batch& b, std::vector<f64>& neg, std::vector<f64>& pos, KernelStats* st) override { std::vector<f64> a, c, f; realized_all(b, a, c, neg, pos, f, st, -1); }
140+ std::vector<f64> sum_fourth(const Batch& b, KernelStats* st) override { std::vector<f64> a, c, d, e, f; realized_all(b, a, c, d, e, f, st, 4); return f; }
141+ /// All five realized accumulators in a single fused pass (used by est::realized_measures_batch when the backend is Metal).
142+ void realized_fused(const Batch& b, std::vector<f64>& s2, std::vector<f64>& adj, std::vector<f64>& neg, std::vector<f64>& pos, std::vector<f64>& s4, KernelStats* st) { realized_all(b, s2, adj, neg, pos, s4, st, -1); }
143+
144+ std::vector<f64> rolling_sum_squares(const Batch& b, usize window, KernelStats* st) override {
145+ Timer t; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
146+ const usize n = b.values.size(), S = b.n_series(); std::vector<f64> out(n, std::numeric_limits<f64>::quiet_NaN()); if (n == 0) return out;
147+ RAID_REQUIRE(window <= 65536, unsupported, "metal rolling window too large for the direct kernel");
148+ @autoreleasepool {
149+ id<MTLBuffer> bx = buffer(n * 4, &allocs), bb = buffer(S * 4, &allocs), bs = buffer(n * 4, &allocs), bo = buffer(n * 4, &allocs);
150+ float* xf = static_cast<float*>(bx.contents); for (usize i = 0; i < n; ++i) xf[i] = static_cast<float>(b.values[i]);
151+ uint32_t* sb = static_cast<uint32_t*>(bb.contents); uint32_t* so = static_cast<uint32_t*>(bs.contents);
152+ for (usize s = 0; s < S; ++s) { sb[s] = static_cast<uint32_t>(b.offsets[s]); for (i64 i = b.offsets[s]; i < b.offsets[s + 1]; ++i) so[i] = static_cast<uint32_t>(s); }
153+ const uint32_t w = static_cast<uint32_t>(window), nn = static_cast<uint32_t>(n);
154+ id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder]; auto p = pipeline("rolling_sumsq");
155+ [enc setComputePipelineState:p]; [enc setBuffer:bx offset:0 atIndex:0]; [enc setBuffer:bb offset:0 atIndex:1]; [enc setBuffer:bs offset:0 atIndex:2]; [enc setBuffer:bo offset:0 atIndex:3]; [enc setBytes:&w length:4 atIndex:4]; [enc setBytes:&nn length:4 atIndex:5];
156+ [enc dispatchThreads:MTLSizeMake(n, 1, 1) threadsPerThreadgroup:MTLSizeMake(std::min<usize>(256, p.maxTotalThreadsPerThreadgroup), 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
157+ RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
158+ const float* o = static_cast<const float*>(bo.contents); for (usize i = 0; i < n; ++i) out[i] = o[i];
159+ finish(st, "rolling_sumsq", t, cb, n * 8, n * 8, allocs, library_source_); release(bx); release(bb); release(bs); release(bo);
160+ }
161+ return out;
162+ }
163+
164+ 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 {
165+ Timer t; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
166+ RAID_REQUIRE(k >= 1 && k <= 8, unsupported, "metal regression_stats supports k <= 8");
167+ usize S; auto chunks = make_chunks(offsets, chunk_size(), &S); const usize rows = y.size(), C = chunks.size(); const usize ntri = k * (k + 1) / 2, nstat = ntri + k + 1;
168+ xtx.assign(S * k * k, 0); xty.assign(S * k, 0); yty.assign(S, 0); if (rows == 0) return;
169+ @autoreleasepool {
170+ id<MTLBuffer> bX = buffer(rows * k * 4, &allocs), bY = buffer(rows * 4, &allocs), bc = buffer(C * sizeof(Chunk), &allocs), bo = buffer(C * nstat * 8, &allocs);
171+ float* xf = static_cast<float*>(bX.contents); for (usize i = 0; i < rows * k; ++i) xf[i] = static_cast<float>(x[i]); float* yf = static_cast<float*>(bY.contents); for (usize i = 0; i < rows; ++i) yf[i] = static_cast<float>(y[i]);
172+ std::memcpy(bc.contents, chunks.data(), C * sizeof(Chunk)); const uint32_t kk = static_cast<uint32_t>(k);
173+ id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder]; auto p = pipeline("regression_stats");
174+ [enc setComputePipelineState:p]; [enc setBuffer:bX offset:0 atIndex:0]; [enc setBuffer:bY offset:0 atIndex:1]; [enc setBuffer:bc offset:0 atIndex:2]; [enc setBuffer:bo offset:0 atIndex:3]; [enc setBytes:&kk length:4 atIndex:4];
175+ [enc dispatchThreadgroups:MTLSizeMake(C, 1, 1) threadsPerThreadgroup:MTLSizeMake(256, 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
176+ RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
177+ const float* o = static_cast<const float*>(bo.contents);
178+ for (usize c = 0; c < C; ++c) { const usize s = chunks[c].seg; usize q = 0;
179+ for (usize a = 0; a < k; ++a) for (usize bcol = a; bcol < k; ++bcol) { const f64 v = static_cast<f64>(o[(c * nstat + q) * 2]) + o[(c * nstat + q) * 2 + 1]; xtx[s * k * k + a * k + bcol] += v; if (a != bcol) xtx[s * k * k + bcol * k + a] += v; ++q; }
180+ for (usize a = 0; a < k; ++a) { xty[s * k + a] += static_cast<f64>(o[(c * nstat + q) * 2]) + o[(c * nstat + q) * 2 + 1]; ++q; }
181+ yty[s] += static_cast<f64>(o[(c * nstat + q) * 2]) + o[(c * nstat + q) * 2 + 1]; }
182+ finish(st, "regression_stats", t, cb, rows * (k + 1) * 8, S * (k * k + k + 1) * 8, allocs, library_source_ + " chunks=" + std::to_string(C)); release(bX); release(bY); release(bc); release(bo);
183+ }
184+ }
185+
186+ std::vector<f64> garch_nll_batch(const Batch& r, std::span<const f64> params, usize n_params, usize p, int model_code, f64 nu_fixed, KernelStats* st) override {
187+ Timer t; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
188+ const usize S = r.n_series(), n = r.values.size(), J = S * n_params; std::vector<f64> out(J, std::numeric_limits<f64>::infinity()); if (J == 0) return out;
189+ RAID_REQUIRE(p == 4 || p == 5 || p == 6, invalid_argument, "garch params per vector must be 4 (garch), 5 (gjr/egarch) or 6 (with nu)");
190+ @autoreleasepool {
191+ id<MTLBuffer> br = buffer(n * 4, &allocs), bb = buffer(S * 4, &allocs), be = buffer(S * 4, &allocs), bp = buffer(J * 8 * 4, &allocs), bv = buffer(S * 4, &allocs), bj = buffer(J * 16, &allocs), bo = buffer(J * 4, &allocs);
192+ float* rf = static_cast<float*>(br.contents); for (usize i = 0; i < n; ++i) rf[i] = static_cast<float>(r.values[i]);
193+ uint32_t* sb = static_cast<uint32_t*>(bb.contents); uint32_t* se = static_cast<uint32_t*>(be.contents); float* v0 = static_cast<float*>(bv.contents);
194+ float* pf = static_cast<float*>(bp.contents); uint32_t* jb = static_cast<uint32_t*>(bj.contents);
195+ const bool gjr = model_code == 1, eg = model_code == 2;
196+ for (usize s = 0; s < S; ++s) { sb[s] = static_cast<uint32_t>(r.offsets[s]); se[s] = static_cast<uint32_t>(r.offsets[s + 1]);
197+ for (usize q = 0; q < n_params; ++q) { const f64* th = params.data() + (s * n_params + q) * p; const usize j = s * n_params + q;
198+ // layout [mu, omega, alpha, gamma, beta, nu]; FP64 initial variance per (series, mu) computed on the CPU → we use mu of the first param set per series
199+ pf[j * 8 + 0] = static_cast<float>(th[0]); pf[j * 8 + 1] = static_cast<float>(th[1]); pf[j * 8 + 2] = static_cast<float>(th[2]);
200+ pf[j * 8 + 3] = (gjr || eg) ? static_cast<float>(th[3]) : 0.0f; pf[j * 8 + 4] = static_cast<float>(th[(gjr || eg) ? 4 : 3]);
201+ const bool has_nu = p == 6 || (p == 5 && !gjr && !eg); const f64 nu = has_nu ? th[p - 1] : (nu_fixed > 2 ? nu_fixed : 0.0); pf[j * 8 + 5] = static_cast<float>(nu);
202+ // FP64 density constants (Metal has no lgamma): kconst and E|z|
203+ const f64 pi = 3.14159265358979323846;
204+ if (nu > 2.0) { const f64 lg = std::lgamma((nu + 1) / 2) - std::lgamma(nu / 2); pf[j * 8 + 6] = static_cast<float>(lg - 0.5 * std::log(pi * (nu - 2.0))); pf[j * 8 + 7] = static_cast<float>(2.0 * std::sqrt(nu - 2.0) * std::exp(lg) / ((nu - 1.0) * std::sqrt(pi))); }
205+ else { pf[j * 8 + 6] = static_cast<float>(-0.5 * std::log(2.0 * pi)); pf[j * 8 + 7] = static_cast<float>(std::sqrt(2.0 / pi)); }
206+ jb[j * 4 + 0] = static_cast<uint32_t>(s); jb[j * 4 + 1] = static_cast<uint32_t>(j); jb[j * 4 + 2] = 0; jb[j * 4 + 3] = 0; }
207+ // initial variance: sample variance of (r − mu) with mu from the first parameter vector of the series (FP64)
208+ const f64 mu = params[(s * n_params) * p]; auto x = r.series(s); f64 m1 = 0; for (f64 v : x) m1 += v - mu; m1 /= static_cast<f64>(x.size()); f64 var0 = 0; for (f64 v : x) var0 += (v - mu - m1) * (v - mu - m1); var0 /= static_cast<f64>(x.size()); v0[s] = static_cast<float>(var0); }
209+ const uint32_t model = static_cast<uint32_t>(model_code), nj = static_cast<uint32_t>(J);
210+ id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder]; auto pl = pipeline("garch_nll");
211+ [enc setComputePipelineState:pl]; [enc setBuffer:br offset:0 atIndex:0]; [enc setBuffer:bb offset:0 atIndex:1]; [enc setBuffer:be offset:0 atIndex:2]; [enc setBuffer:bp offset:0 atIndex:3]; [enc setBuffer:bv offset:0 atIndex:4]; [enc setBuffer:bj offset:0 atIndex:5]; [enc setBuffer:bo offset:0 atIndex:6]; [enc setBytes:&model length:4 atIndex:7]; [enc setBytes:&nj length:4 atIndex:8];
212+ [enc dispatchThreads:MTLSizeMake(J, 1, 1) threadsPerThreadgroup:MTLSizeMake(std::min<usize>(64, pl.maxTotalThreadsPerThreadgroup), 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
213+ RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
214+ const float* o = static_cast<const float*>(bo.contents); for (usize j = 0; j < J; ++j) out[j] = o[j];
215+ finish(st, "garch_nll", t, cb, n * 8 + J * p * 8, J * 8, allocs, library_source_ + " jobs=" + std::to_string(J)); release(br); release(bb); release(be); release(bp); release(bv); release(bj); release(bo);
216+ }
217+ return out;
218+ }
219+ f64 compile_seconds() const { return compile_seconds_; }
220+ std::string library_source() const { return library_source_; }
221+
222+ private:
223+ PrecisionPolicy policy_; id<MTLDevice> device_ = nil; id<MTLCommandQueue> queue_ = nil; id<MTLLibrary> library_ = nil; bool available_{false};
224+ std::map<std::string, id<MTLComputePipelineState>> pipelines_; std::map<usize, std::vector<id<MTLBuffer>>> pool_; std::mutex m_; usize total_allocs_{0};
225+ std::string library_source_{"none"}; f64 compile_seconds_{0};
226+};
227+
228+std::shared_ptr<Backend> make_metal_backend(const PrecisionPolicy& policy) {
229+ static std::mutex m; static std::weak_ptr<MetalBackend> cached; std::lock_guard<std::mutex> lk(m);
230+ if (auto p = cached.lock()) return p; // one device/queue/library per process (pipelines and buffers are reused)
231+ auto p = std::make_shared<MetalBackend>(policy); cached = p; if (!p->available()) return nullptr; return p;
232+}
233+
234+/// Fused realized accumulators when the backend is Metal (falls back to the four generic calls otherwise).
235+bool realized_fused(Backend& be, const Batch& b, std::vector<f64>& s2, std::vector<f64>& adj, std::vector<f64>& neg, std::vector<f64>& pos, std::vector<f64>& s4, KernelStats* st) {
236+ if (auto* m = dynamic_cast<MetalBackend*>(&be)) { m->realized_fused(b, s2, adj, neg, pos, s4, st); return true; }
237+ return false;
238+}
239+f64 metal_compile_seconds(Backend& be) { if (auto* m = dynamic_cast<MetalBackend*>(&be)) return m->compile_seconds(); return 0; }
240+
241+} // namespace raid::backend
added backends/metal/shaders/raid_kernels.metal +129 −0
@@ -0,0 +1,129 @@
1+// RAID Metal kernels — FP32 with double-float ("compensated") accumulation.
2+// Apple GPUs have no FP64. Every accumulator is a pair (hi, lo) maintained with error-free transforms (TwoSum /
3+// Fast2Sum); partial results per (series, chunk) are written as float2 and combined in FP64 on the CPU in a fixed
4+// order, which makes results deterministic for a given chunk size and threadgroup width.
5+#include <metal_stdlib>
6+using namespace metal;
7+
8+// ---- double-float arithmetic -----------------------------------------------------------------------------------
9+struct df { float hi; float lo; };
10+inline df df_zero() { return df{0.0f, 0.0f}; }
11+inline df two_sum(float a, float b) { float s = a + b; float bb = s - a; float e = (a - (s - bb)) + (b - bb); return df{s, e}; }
12+inline df df_add_f(df a, float b) { df s = two_sum(a.hi, b); s.lo += a.lo; float h = s.hi + s.lo; float l = s.lo - (h - s.hi); return df{h, l}; }
13+inline df df_add(df a, df b) { df s = two_sum(a.hi, b.hi); s.lo += a.lo + b.lo; float h = s.hi + s.lo; float l = s.lo - (h - s.hi); return df{h, l}; }
14+inline float2 df_pack(df a) { return float2(a.hi, a.lo); }
15+
16+// Threadgroup tree reduction of double-float values (deterministic order).
17+inline df tg_reduce(df v, threadgroup df* sh, uint lid, uint n) {
18+ sh[lid] = v; threadgroup_barrier(mem_flags::mem_threadgroup);
19+ for (uint s = n / 2; s > 0; s >>= 1) { if (lid < s) sh[lid] = df_add(sh[lid], sh[lid + s]); threadgroup_barrier(mem_flags::mem_threadgroup); }
20+ return sh[0];
21+}
22+
23+// ---- 1. fused realized-measure accumulators ---------------------------------------------------------------------
24+// One threadgroup per (segment, chunk). Inputs: x (float), offsets (int64 as two uints packed → we use uint offsets,
25+// the host guarantees < 2^32 elements), chunk descriptors (segment id, begin, end). Outputs per chunk: 5 float2:
26+// Σx², Σ|x_i||x_{i-1}| (pairs with i-1 inside the chunk; the host adds the boundary products), Σx²1{x<0}, Σx²1{x>0}, Σx⁴.
27+struct Chunk { uint seg; uint begin; uint end; uint pad; };
28+kernel void realized_fused(device const float* x [[buffer(0)]], device const Chunk* chunks [[buffer(1)]], device float2* out [[buffer(2)]],
29+ uint gid [[threadgroup_position_in_grid]], uint lid [[thread_position_in_threadgroup]], uint tsz [[threads_per_threadgroup]]) {
30+ threadgroup df sh[256];
31+ Chunk c = chunks[gid];
32+ df s2 = df_zero(), sadj = df_zero(), sneg = df_zero(), spos = df_zero(), s4 = df_zero();
33+ for (uint i = c.begin + lid; i < c.end; i += tsz) {
34+ float v = x[i]; float q = v * v;
35+ s2 = df_add_f(s2, q); s4 = df_add_f(s4, q * q);
36+ if (v < 0.0f) sneg = df_add_f(sneg, q); else if (v > 0.0f) spos = df_add_f(spos, q);
37+ if (i > c.begin) sadj = df_add_f(sadj, fabs(v) * fabs(x[i - 1]));
38+ }
39+ df r2 = tg_reduce(s2, sh, lid, tsz); threadgroup_barrier(mem_flags::mem_threadgroup);
40+ df radj = tg_reduce(sadj, sh, lid, tsz); threadgroup_barrier(mem_flags::mem_threadgroup);
41+ df rneg = tg_reduce(sneg, sh, lid, tsz); threadgroup_barrier(mem_flags::mem_threadgroup);
42+ df rpos = tg_reduce(spos, sh, lid, tsz); threadgroup_barrier(mem_flags::mem_threadgroup);
43+ df r4 = tg_reduce(s4, sh, lid, tsz);
44+ if (lid == 0) { out[gid * 5 + 0] = df_pack(r2); out[gid * 5 + 1] = df_pack(radj); out[gid * 5 + 2] = df_pack(rneg); out[gid * 5 + 3] = df_pack(rpos); out[gid * 5 + 4] = df_pack(r4); }
45+}
46+
47+// Single-statistic variants (for the fusion ablation): op 0 = Σx², 1 = Σ|x||x_{-1}|, 2 = Σx² 1{x<0}, 3 = Σx² 1{x>0}, 4 = Σx⁴
48+kernel void segment_reduce(device const float* x [[buffer(0)]], device const Chunk* chunks [[buffer(1)]], device float2* out [[buffer(2)]], constant uint& op [[buffer(3)]],
49+ uint gid [[threadgroup_position_in_grid]], uint lid [[thread_position_in_threadgroup]], uint tsz [[threads_per_threadgroup]]) {
50+ threadgroup df sh[256];
51+ Chunk c = chunks[gid]; df acc = df_zero();
52+ for (uint i = c.begin + lid; i < c.end; i += tsz) {
53+ float v = x[i]; float q = v * v; float t = 0.0f;
54+ switch (op) { case 0: t = q; break; case 1: t = (i > c.begin) ? fabs(v) * fabs(x[i - 1]) : 0.0f; break; case 2: t = v < 0.0f ? q : 0.0f; break; case 3: t = v > 0.0f ? q : 0.0f; break; default: t = q * q; }
55+ acc = df_add_f(acc, t);
56+ }
57+ df r = tg_reduce(acc, sh, lid, tsz);
58+ if (lid == 0) out[gid] = df_pack(r);
59+}
60+
61+// ---- 2. rolling sum of squares (direct window, O(w) per output) ---------------------------------------------------
62+// out[i] = Σ_{j=i-w+1..i} x_j² within the segment; NaN where the window is incomplete. One thread per element.
63+kernel void rolling_sumsq(device const float* x [[buffer(0)]], device const uint* seg_begin [[buffer(1)]], device const uint* seg_of [[buffer(2)]], device float* out [[buffer(3)]],
64+ constant uint& w [[buffer(4)]], constant uint& n [[buffer(5)]], uint i [[thread_position_in_grid]]) {
65+ if (i >= n) return;
66+ uint b = seg_begin[seg_of[i]];
67+ if (i + 1 < b + w) { out[i] = NAN; return; }
68+ df acc = df_zero();
69+ for (uint j = i + 1 - w; j <= i; ++j) { float v = x[j]; acc = df_add_f(acc, v * v); }
70+ out[i] = acc.hi + acc.lo;
71+}
72+
73+// ---- 3. regression sufficient statistics (X'X upper triangle, X'y, y'y) per chunk ------------------------------
74+// X row-major (rows × k), k ≤ 8. One threadgroup per chunk; each thread accumulates its rows in double-float registers,
75+// then a deterministic tree reduction per statistic.
76+#define RAID_MAXK 8
77+kernel void regression_stats(device const float* X [[buffer(0)]], device const float* y [[buffer(1)]], device const Chunk* chunks [[buffer(2)]], device float2* out [[buffer(3)]],
78+ constant uint& k [[buffer(4)]], uint gid [[threadgroup_position_in_grid]], uint lid [[thread_position_in_threadgroup]], uint tsz [[threads_per_threadgroup]]) {
79+ threadgroup df sh[256];
80+ Chunk c = chunks[gid];
81+ const uint ntri = k * (k + 1) / 2; const uint nstat = ntri + k + 1;
82+ df acc[RAID_MAXK * (RAID_MAXK + 1) / 2 + RAID_MAXK + 1];
83+ for (uint s = 0; s < nstat; ++s) acc[s] = df_zero();
84+ for (uint r = c.begin + lid; r < c.end; r += tsz) {
85+ float xr[RAID_MAXK]; for (uint a = 0; a < k; ++a) xr[a] = X[r * k + a];
86+ float yr = y[r]; uint s = 0;
87+ for (uint a = 0; a < k; ++a) for (uint b = a; b < k; ++b) { acc[s] = df_add_f(acc[s], xr[a] * xr[b]); ++s; }
88+ for (uint a = 0; a < k; ++a) { acc[s] = df_add_f(acc[s], xr[a] * yr); ++s; }
89+ acc[s] = df_add_f(acc[s], yr * yr);
90+ }
91+ for (uint s = 0; s < nstat; ++s) { df r = tg_reduce(acc[s], sh, lid, tsz); if (lid == 0) out[gid * nstat + s] = df_pack(r); threadgroup_barrier(mem_flags::mem_threadgroup); }
92+}
93+
94+// ---- 4. batched GARCH-family negative log-likelihood ---------------------------------------------------------------
95+// One thread per (series, parameter vector). params layout per thread (stride 8): [mu, omega, alpha, gamma, beta, nu,
96+// kconst, eabs] — gamma=0 for GARCH, nu<=2 → Gaussian; kconst = density constant and eabs = E|z| are computed in FP64
97+// on the host (Metal has no lgamma). var0 = FP64-computed initial variance passed from the host (mixed precision).
98+// model: 0 garch, 1 gjr, 2 egarch. Returns NLL (double-float summed) as float; +INF for invalid parameters.
99+struct GarchJob { uint seg; uint pidx; uint pad0; uint pad1; };
100+kernel void garch_nll(device const float* r [[buffer(0)]], device const uint* seg_begin [[buffer(1)]], device const uint* seg_end [[buffer(2)]], device const float* params [[buffer(3)]],
101+ device const float* var0 [[buffer(4)]], device const GarchJob* jobs [[buffer(5)]], device float* out [[buffer(6)]], constant uint& model [[buffer(7)]], constant uint& njobs [[buffer(8)]],
102+ uint t [[thread_position_in_grid]]) {
103+ if (t >= njobs) return;
104+ GarchJob j = jobs[t]; const uint b = seg_begin[j.seg], e = seg_end[j.seg]; const uint n = e - b;
105+ const float mu = params[j.pidx * 8 + 0], omega = params[j.pidx * 8 + 1], alpha = params[j.pidx * 8 + 2], gamma = params[j.pidx * 8 + 3], beta = params[j.pidx * 8 + 4], nu = params[j.pidx * 8 + 5];
106+ const float kconst = params[j.pidx * 8 + 6], Eabs = params[j.pidx * 8 + 7];
107+ const bool student = nu > 2.0f;
108+ bool valid = true;
109+ if (model == 0u) valid = omega > 0.0f && alpha >= 0.0f && beta >= 0.0f && alpha + beta < 1.0f;
110+ else if (model == 1u) valid = omega > 0.0f && alpha >= 0.0f && beta >= 0.0f && alpha + gamma >= 0.0f && alpha + beta + 0.5f * gamma < 1.0f;
111+ else valid = fabs(beta) < 1.0f;
112+ if (!valid || n < 2) { out[t] = INFINITY; return; }
113+ float s2 = var0[j.seg];
114+ df nll = df_zero();
115+ float e_prev = r[b] - mu;
116+ for (uint i = 0; i < n; ++i) {
117+ if (i > 0) {
118+ if (model == 2u) { float z = e_prev / sqrt(s2); float ls = omega + alpha * (fabs(z) - Eabs) + gamma * z + beta * log(s2); s2 = exp(ls); }
119+ else { float e2 = e_prev * e_prev; s2 = omega + (alpha + (e_prev < 0.0f ? gamma : 0.0f)) * e2 + beta * s2; }
120+ if (!(s2 > 0.0f) || !isfinite(s2)) { out[t] = INFINITY; return; }
121+ }
122+ float ev = r[b + i] - mu; float ll;
123+ if (student) ll = kconst - 0.5f * log(s2) - 0.5f * (nu + 1.0f) * log(1.0f + ev * ev / ((nu - 2.0f) * s2)); // no log1p in MSL; argument is O(1)
124+ else ll = kconst - 0.5f * log(s2) - 0.5f * ev * ev / s2;
125+ nll = df_add_f(nll, -ll);
126+ e_prev = ev;
127+ }
128+ out[t] = nll.hi + nll.lo;
129+}
added benchmarks/micro/bench_primitives.cpp +64 −0
@@ -0,0 +1,64 @@
1+// Microbenchmark: batched realized accumulators (fused), rolling sum of squares, regression sufficient statistics
2+// and batched GARCH likelihoods — RAID CPU 1 thread vs CPU all threads vs Metal — over series counts and lengths.
3+// Protocol (charter §8): warm-ups, randomized order, R repetitions, medians + bootstrap CIs computed later from the raw
4+// records; cold-vs-warm shader compilation measured separately; records written as JSON lines with the manifest.
5+#include <cstdio>
6+#include <algorithm>
7+#include <random>
8+#include <fstream>
9+#include <chrono>
10+#include "raid/backend/backend.hpp"
11+#include "raid/core/rng.hpp"
12+#include "raid/core/manifest.hpp"
13+#include "raid/core/thread_pool.hpp"
14+#include "raid/estimators/realized.hpp"
15+using namespace raid; using json = nlohmann::json;
16+
17+namespace {
18+struct Case { std::string kernel; usize S, n; std::string backend; unsigned threads; };
19+f64 now() { return std::chrono::duration<f64>(std::chrono::steady_clock::now().time_since_epoch()).count(); }
20+}
21+
22+int main(int argc, char** argv) {
23+ std::string out = argc > 1 ? argv[1] : "experiments/results/raw/bench_primitives.jsonl"; int reps = argc > 2 ? std::atoi(argv[2]) : 7; const u64 seed = 20261102;
24+ std::vector<usize> Ss = {1, 10, 100, 1000, 10000}, ns = {390, 3900, 39000, 390000};
25+ const usize max_elems = 200'000'000; // memory bound (1.6 GB FP64 + 0.8 GB FP32)
26+ std::vector<Case> cases;
27+ for (auto S : Ss) for (auto n : ns) { if (S * n > max_elems) continue; for (auto k : {"realized_fused", "rolling_sumsq", "regression_stats"}) { cases.push_back({k, S, n, "cpu", 1}); cases.push_back({k, S, n, "cpu", 0}); cases.push_back({k, S, n, "metal", 0}); } }
28+ for (auto S : {1ul, 10ul, 100ul, 1000ul, 10000ul}) for (auto n : {390ul, 3900ul}) for (auto P : {1ul, 64ul}) { if (S * P > 200000) continue; cases.push_back({"garch_nll_P" + std::to_string(P), S, n, "cpu", 1}); cases.push_back({"garch_nll_P" + std::to_string(P), S, n, "cpu", 0}); cases.push_back({"garch_nll_P" + std::to_string(P), S, n, "metal", 0}); }
29+ std::mt19937_64 order(seed); std::shuffle(cases.begin(), cases.end(), order);
30+ ExperimentManifest man; man.experiment = "bench_primitives"; man.created_at_utc = now_utc_iso(); man.raid_version = RAID_VERSION; man.host = host_info(); man.seed = seed; man.parameters = {{"reps", std::to_string(reps)}, {"order", "shuffled"}, {"chunk_max", "4096"}};
31+ std::ofstream f(out, std::ios::app); if (!f.good()) { std::fprintf(stderr, "cannot write %s\n", out.c_str()); return 1; }
32+ auto cpu1 = backend::make_cpu_backend(1), cpuN = backend::make_cpu_backend(0);
33+ std::shared_ptr<backend::Backend> gpu;
34+#ifdef RAID_ENABLE_METAL
35+ { auto t0 = now(); gpu = backend::make_metal_backend(); if (gpu) { backend::KernelStats st; std::vector<f64> v(390, 1e-3); std::vector<i64> off = {0, 390}; backend::Batch b{v, off}; std::vector<f64> a, b2, c, d, e; backend::realized_fused(*gpu, b, a, b2, c, d, e, &st);
36+ json rec = {{"kernel", "shader_library_load"}, {"backend", "metal"}, {"compiled_now", st.compiled_now}, {"seconds", backend::metal_compile_seconds(*gpu)}, {"wall_first_call", st.wall_seconds}, {"note", st.note}, {"t", now() - t0}}; f << rec.dump() << "\n"; } }
37+#endif
38+ usize done = 0;
39+ for (auto& c : cases) {
40+ if (c.backend == "metal" && !gpu) continue;
41+ auto be = c.backend == "metal" ? gpu : (c.threads == 1 ? cpu1 : cpuN);
42+ Rng r(seed, c.S * 1000003 + c.n); std::vector<f64> v(c.S * c.n); for (auto& x : v) x = 1e-3 * r.student_t(5.0); std::vector<i64> off(c.S + 1); for (usize s = 0; s <= c.S; ++s) off[s] = static_cast<i64>(s * c.n);
43+ backend::Batch b{v, off};
44+ std::vector<f64> X, y; usize k = 4; usize P = 1; std::vector<f64> params;
45+ if (c.kernel == "regression_stats") { X.resize(c.S * c.n * k); y.resize(c.S * c.n); for (usize i = 0; i < c.S * c.n; ++i) { X[i * k] = 1; for (usize j = 1; j < k; ++j) X[i * k + j] = r.normal(); y[i] = r.normal(); } }
46+ if (c.kernel.rfind("garch_nll_P", 0) == 0) { P = std::stoul(c.kernel.substr(11)); params.resize(c.S * P * 4); for (usize s = 0; s < c.S; ++s) for (usize q = 0; q < P; ++q) { f64* th = params.data() + (s * P + q) * 4; th[0] = 0; th[1] = 1e-7; th[2] = 0.05 + 0.001 * q; th[3] = 0.9; } }
47+ auto run = [&](backend::KernelStats& st) {
48+ if (c.kernel == "realized_fused") { std::vector<f64> a, b2, cc, d, e;
49+#ifdef RAID_ENABLE_METAL
50+ if (c.backend == "metal") { backend::realized_fused(*be, b, a, b2, cc, d, e, &st); return; }
51+#endif
52+ backend::KernelStats s1, s2, s3, s4; a = be->sum_squares(b, &s1); b2 = be->sum_abs_adjacent(b, &s2); be->semivariances(b, cc, d, &s3); e = be->sum_fourth(b, &s4); st = s1; st.wall_seconds = s1.wall_seconds + s2.wall_seconds + s3.wall_seconds + s4.wall_seconds; }
53+ else if (c.kernel == "rolling_sumsq") { auto o = be->rolling_sum_squares(b, 22, &st); }
54+ else if (c.kernel == "regression_stats") { std::vector<f64> a, bb, cc; be->regression_stats(X, y, off, k, a, bb, cc, &st); }
55+ else { auto o = be->garch_nll_batch(b, params, P, 4, 0, 0.0, &st); }
56+ };
57+ backend::KernelStats warm; for (int w = 0; w < 2; ++w) run(warm);
58+ for (int rep = 0; rep < reps; ++rep) { backend::KernelStats st; const u64 rss0 = peak_rss_bytes(); const f64 t0 = now(); run(st); const f64 wall = now() - t0;
59+ json rec = {{"kernel", c.kernel}, {"S", c.S}, {"n", c.n}, {"elements", c.S * c.n}, {"backend", c.backend}, {"threads", c.backend == "metal" ? 0u : (c.threads == 1 ? 1u : ThreadPool::global().size())}, {"rep", rep}, {"wall_seconds", wall}, {"kernel_wall_seconds", st.wall_seconds}, {"device_seconds", st.device_seconds}, {"bytes_in", st.bytes_in}, {"buffer_allocations", st.buffer_allocations}, {"precision", std::string(to_string(st.precision))}, {"note", st.note}, {"peak_rss_before", rss0}, {"peak_rss_after", peak_rss_bytes()}, {"P", P}, {"k", k}, {"window", 22}};
60+ f << rec.dump() << "\n"; }
61+ f.flush(); ++done; std::fprintf(stderr, "[%zu/%zu] %s S=%zu n=%zu %s t=%u\n", done, cases.size(), c.kernel.c_str(), c.S, c.n, c.backend.c_str(), c.threads);
62+ }
63+ man.save(out + ".manifest.json"); std::fprintf(stderr, "wrote %s (+ manifest)\n", out.c_str()); return 0;
64+}
added experiments/scripts/analyze_bench_primitives.py +71 −0
@@ -0,0 +1,71 @@
1+#!/usr/bin/env python3
2+"""Break-even analysis of bench_primitives.jsonl: medians, bootstrap 95 % CIs, speedups T_ref/T_metal, break-even
3+elements per kernel, CPU-wins table; writes experiments/results/bench_primitives_summary.{json,csv} and figures.
4+Usage: analyze_bench_primitives.py [raw.jsonl] [out_dir]
5+"""
6+import json, os, sys, collections
7+import numpy as np
8+
9+raw = sys.argv[1] if len(sys.argv) > 1 else "experiments/results/raw/bench_primitives.jsonl"
10+out_dir = sys.argv[2] if len(sys.argv) > 2 else "experiments/results"
11+recs = [json.loads(l) for l in open(raw) if l.strip()]
12+lib = [r for r in recs if r.get("kernel") == "shader_library_load"]
13+recs = [r for r in recs if r.get("kernel") != "shader_library_load"]
14+rng = np.random.default_rng(20261102)
15+
16+def boot_ci(x, n=2000):
17+ x = np.asarray(x); meds = [np.median(rng.choice(x, x.size, replace=True)) for _ in range(n)]
18+ return float(np.median(x)), float(np.percentile(meds, 2.5)), float(np.percentile(meds, 97.5))
19+
20+groups = collections.defaultdict(list)
21+for r in recs:
22+ key = (r["kernel"], r["S"], r["n"], r["backend"], r["threads"])
23+ groups[key].append(r)
24+rows = []
25+for (kernel, S, n, backend, threads), rs in sorted(groups.items()):
26+ w = [r["wall_seconds"] for r in rs]; med, lo, hi = boot_ci(w)
27+ dev = float(np.median([r["device_seconds"] for r in rs])) if backend == "metal" else 0.0
28+ rows.append({"kernel": kernel, "S": S, "n": n, "elements": S * n, "backend": backend, "threads": threads, "median_s": med, "ci_lo": lo, "ci_hi": hi,
29+ "device_s": dev, "elements_per_s": S * n / med, "reps": len(rs), "precision": rs[0]["precision"], "peak_rss_mb": max(r["peak_rss_after"] for r in rs) / 1e6})
30+# speedups vs CPU multi-thread (reference) and CPU single-thread
31+by = {(r["kernel"], r["S"], r["n"], r["backend"], r["threads"] == 1): r for r in rows}
32+summary = []
33+for r in rows:
34+ if r["backend"] != "metal": continue
35+ cpuN = by.get((r["kernel"], r["S"], r["n"], "cpu", False)); cpu1 = by.get((r["kernel"], r["S"], r["n"], "cpu", True))
36+ summary.append({**r, "speedup_vs_cpu_mt": cpuN["median_s"] / r["median_s"] if cpuN else None, "speedup_vs_cpu_1t": cpu1["median_s"] / r["median_s"] if cpu1 else None,
37+ "cpu_mt_median_s": cpuN["median_s"] if cpuN else None, "cpu_1t_median_s": cpu1["median_s"] if cpu1 else None,
38+ "metal_wins_vs_mt": bool(cpuN and r["ci_hi"] < cpuN["ci_lo"]), "cpu_wins_vs_mt": bool(cpuN and cpuN["ci_hi"] < r["ci_lo"])})
39+# break-even: smallest element count where Metal beats CPU-MT (CI-separated) per kernel
40+breakeven = {}
41+for k in sorted(set(r["kernel"] for r in summary)):
42+ wins = sorted([r["elements"] for r in summary if r["kernel"] == k and r["metal_wins_vs_mt"]])
43+ losses = sorted([r["elements"] for r in summary if r["kernel"] == k and r["cpu_wins_vs_mt"]])
44+ breakeven[k] = {"smallest_metal_win_elements": wins[0] if wins else None, "largest_cpu_win_elements": losses[-1] if losses else None,
45+ "max_speedup_vs_mt": max((r["speedup_vs_cpu_mt"] or 0) for r in summary if r["kernel"] == k), "n_cases": sum(1 for r in summary if r["kernel"] == k)}
46+os.makedirs(out_dir, exist_ok=True)
47+json.dump({"source": raw, "rows": rows, "metal_summary": summary, "breakeven": breakeven, "shader_library": lib}, open(os.path.join(out_dir, "bench_primitives_summary.json"), "w"), indent=1)
48+import csv
49+with open(os.path.join(out_dir, "bench_primitives_summary.csv"), "w", newline="") as f:
50+ w = csv.DictWriter(f, fieldnames=list(rows[0].keys())); w.writeheader(); w.writerows(rows)
51+print("kernel S n elements cpu1(ms) cpuMT(ms) metal(ms) dev(ms) speedup_MT verdict")
52+for r in sorted(summary, key=lambda r: (r["kernel"], r["elements"], r["S"])):
53+ v = "METAL" if r["metal_wins_vs_mt"] else ("CPU" if r["cpu_wins_vs_mt"] else "tie")
54+ print(f"{r['kernel']:20s} {r['S']:6d} {r['n']:7d} {r['elements']:11d} {1e3*(r['cpu_1t_median_s'] or 0):9.3f} {1e3*(r['cpu_mt_median_s'] or 0):9.3f} {1e3*r['median_s']:9.3f} {1e3*r['device_s']:8.3f} {r['speedup_vs_cpu_mt'] or 0:10.2f} {v}")
55+print("break-even:", json.dumps(breakeven, indent=1))
56+if lib: print("shader library:", lib)
57+try:
58+ import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
59+ kernels = sorted(set(r["kernel"] for r in rows)); fig, axes = plt.subplots(1, len(kernels), figsize=(4.2 * len(kernels), 3.6), squeeze=False)
60+ for ax, k in zip(axes[0], kernels):
61+ for backend, thr, lab, col in [("cpu", True, "CPU 1 thread", "#9aa0a6"), ("cpu", False, "CPU multicore", "#1f4e79"), ("metal", False, "Metal (FP32 df)", "#c0392b")]:
62+ pts = sorted([(r["elements"], r["elements_per_s"], r["ci_lo"], r["ci_hi"], r["median_s"]) for r in rows if r["kernel"] == k and r["backend"] == backend and (r["threads"] == 1) == thr])
63+ if not pts: continue
64+ x = [p[0] for p in pts]; y = [p[1] for p in pts]
65+ ax.plot(x, y, "o-", ms=3, lw=1.2, label=lab, color=col)
66+ ax.set_xscale("log"); ax.set_yscale("log"); ax.set_title(k, fontsize=9); ax.set_xlabel("elements (S×n)"); ax.grid(alpha=0.3)
67+ axes[0][0].set_ylabel("elements / s (median)"); axes[0][0].legend(fontsize=7)
68+ fig.tight_layout(); fig.savefig(os.path.join(out_dir, "fig_bench_primitives_throughput.pdf")); fig.savefig(os.path.join(out_dir, "fig_bench_primitives_throughput.png"), dpi=150)
69+ print("figure written")
70+except Exception as e:
71+ print("figure skipped:", e)
modified include/raid/backend/backend.hpp +9 −0
@@ -69,3 +69,12 @@ std::shared_ptr<Backend> resolve(const ExecutionPlan& plan, std::string* reason
69 std::vector<Device> list_devices(); 69 std::vector<Device> list_devices();
70 70
71 } // namespace raid::backend 71 } // namespace raid::backend
72+
73+namespace raid::backend {
74+#ifdef RAID_ENABLE_METAL
75+/// Metal-only fused realized kernel (one pass for Σx², Σ|x||x₋₁|, Σx²1{x<0}, Σx²1{x>0}, Σx⁴). Returns false if `be` is not Metal.
76+bool realized_fused(Backend& be, const Batch& b, std::vector<f64>& s2, std::vector<f64>& adj, std::vector<f64>& neg, std::vector<f64>& pos, std::vector<f64>& s4, KernelStats* st = nullptr);
77+/// Seconds spent compiling/loading the shader library (0 if not Metal or not yet loaded).
78+f64 metal_compile_seconds(Backend& be);
79+#endif
80+}
modified include/raid/estimators/realized.hpp +1 −0
@@ -24,6 +24,7 @@ namespace raid::est {
24 24
25 struct RealizedOptions { 25 struct RealizedOptions {
26 bool bpv_finite_sample{true}; 26 bool bpv_finite_sample{true};
27+ bool fuse_kernels{true}; // Metal: one fused pass for all accumulators (ablation: false → 4 kernels)
27 bool annualize{false}; // multiply variances by `periods_per_year` 28 bool annualize{false}; // multiply variances by `periods_per_year`
28 f64 periods_per_year{252.0}; 29 f64 periods_per_year{252.0};
29 ExecutionPlan plan{}; 30 ExecutionPlan plan{};
modified src/estimators/realized.cpp +12 −7
@@ -11,13 +11,18 @@ f64 mu_p(f64 p) { return std::pow(2.0, p / 2.0) * std::tgamma((p + 1.0) / 2.0) /
11 RealizedMeasures realized_measures_batch(const backend::Batch& b, const RealizedOptions& opt) { 11 RealizedMeasures realized_measures_batch(const backend::Batch& b, const RealizedOptions& opt) {
12 RealizedMeasures m; const usize S = b.n_series(); 12 RealizedMeasures m; const usize S = b.n_series();
13 std::string reason; auto be = backend::resolve(opt.plan, &reason); 13 std::string reason; auto be = backend::resolve(opt.plan, &reason);
14− backend::KernelStats st1, st2, st3, st4; 14+ backend::KernelStats st1, st2, st3, st4; std::vector<f64> adj, s4; bool fused = false;
15− m.rv = be->sum_squares(b, &st1); 15+#ifdef RAID_ENABLE_METAL
16− auto adj = be->sum_abs_adjacent(b, &st2); 16+ if (opt.fuse_kernels) fused = backend::realized_fused(*be, b, m.rv, adj, m.rsv_neg, m.rsv_pos, s4, &st1);
17− be->semivariances(b, m.rsv_neg, m.rsv_pos, &st3); 17+#endif
18− auto s4 = be->sum_fourth(b, &st4); 18+ if (!fused) {
19− m.stats = st1; m.stats.kernel = "realized_measures"; m.stats.wall_seconds = st1.wall_seconds + st2.wall_seconds + st3.wall_seconds + st4.wall_seconds; 19+ m.rv = be->sum_squares(b, &st1);
20− m.stats.device_seconds = st1.device_seconds + st2.device_seconds + st3.device_seconds + st4.device_seconds; m.stats.note = reason; 20+ adj = be->sum_abs_adjacent(b, &st2);
21+ be->semivariances(b, m.rsv_neg, m.rsv_pos, &st3);
22+ s4 = be->sum_fourth(b, &st4);
23+ }
24+ m.stats = st1; m.stats.kernel = fused ? "realized_fused" : "realized_measures(4 kernels)"; m.stats.wall_seconds = st1.wall_seconds + st2.wall_seconds + st3.wall_seconds + st4.wall_seconds;
25+ m.stats.device_seconds = st1.device_seconds + st2.device_seconds + st3.device_seconds + st4.device_seconds; m.stats.note = reason + (m.stats.note.empty() ? "" : " | " + m.stats.note);
21 m.bpv.resize(S); m.rq.resize(S); m.jump.resize(S); m.tripv.resize(S); m.n_returns.resize(S); m.daily_return.resize(S); 26 m.bpv.resize(S); m.rq.resize(S); m.jump.resize(S); m.tripv.resize(S); m.n_returns.resize(S); m.daily_return.resize(S);
22 const f64 mu43 = mu_p(4.0 / 3.0); 27 const f64 mu43 = mu_p(4.0 / 3.0);
23 for (usize s = 0; s < S; ++s) { 28 for (usize s = 0; s < S; ++s) {
added tests/unit/test_metal_backend.cpp +79 −0
@@ -0,0 +1,79 @@
1+// Metal backend: accuracy vs the FP64 CPU reference, determinism, chunking, fallback. Skipped when Metal is absent.
2+#include <catch2/catch_test_macros.hpp>
3+#include <catch2/matchers/catch_matchers_floating_point.hpp>
4+#include "raid/backend/backend.hpp"
5+#include "raid/core/rng.hpp"
6+#include "raid/estimators/realized.hpp"
7+#include "raid/estimators/ols.hpp"
8+#include "raid/models/garch.hpp"
9+#include <cmath>
10+using namespace raid; using Catch::Matchers::WithinRel; using Catch::Matchers::WithinAbs;
11+
12+#ifdef RAID_ENABLE_METAL
13+static std::shared_ptr<backend::Backend> metal() { auto m = backend::make_metal_backend(); return (m && m->available()) ? m : nullptr; }
14+
15+static void make_batch(usize S, usize n, std::vector<f64>& v, std::vector<i64>& off, u64 seed = 1) {
16+ Rng r(seed); v.resize(S * n); off.resize(S + 1); for (usize s = 0; s <= S; ++s) off[s] = static_cast<i64>(s * n);
17+ for (auto& x : v) x = 1e-3 * r.student_t(5.0);
18+}
19+
20+TEST_CASE("Metal realized accumulators match FP64 CPU within 1e-6 (fused and single kernels, chunked)", "[metal]") {
21+ auto m = metal(); if (!m) { SKIP("no Metal device"); }
22+ auto cpu = backend::make_cpu_backend();
23+ for (usize n : {390ul, 3900ul, 39000ul}) {
24+ std::vector<f64> v; std::vector<i64> off; make_batch(7, n, v, off, n); backend::Batch b{v, off};
25+ backend::KernelStats st; std::vector<f64> s2, adj, neg, pos, s4; REQUIRE(backend::realized_fused(*m, b, s2, adj, neg, pos, s4, &st));
26+ auto c2 = cpu->sum_squares(b), cadj = cpu->sum_abs_adjacent(b), c4 = cpu->sum_fourth(b); std::vector<f64> cn, cp; cpu->semivariances(b, cn, cp);
27+ for (usize s = 0; s < 7; ++s) { INFO("n=" << n << " s=" << s); CHECK_THAT(s2[s], WithinRel(c2[s], 1e-6)); CHECK_THAT(adj[s], WithinRel(cadj[s], 1e-6)); CHECK_THAT(neg[s], WithinRel(cn[s], 1e-6)); CHECK_THAT(pos[s], WithinRel(cp[s], 1e-6)); CHECK_THAT(s4[s], WithinRel(c4[s], 1e-5)); }
28+ CHECK(st.backend == Backend::metal); CHECK(st.precision == Precision::fp32); CHECK(st.device_seconds > 0);
29+ auto g2 = m->sum_squares(b); auto gadj = m->sum_abs_adjacent(b); for (usize s = 0; s < 7; ++s) { CHECK_THAT(g2[s], WithinRel(c2[s], 1e-6)); CHECK_THAT(gadj[s], WithinRel(cadj[s], 1e-6)); }
30+ }
31+}
32+TEST_CASE("Metal results are bitwise deterministic across runs", "[metal]") {
33+ auto m = metal(); if (!m) { SKIP("no Metal device"); }
34+ std::vector<f64> v; std::vector<i64> off; make_batch(50, 3900, v, off, 9); backend::Batch b{v, off};
35+ std::vector<f64> a2, aa, an, ap, a4, b2, ba, bn, bp, b4; backend::realized_fused(*m, b, a2, aa, an, ap, a4); backend::realized_fused(*m, b, b2, ba, bn, bp, b4);
36+ CHECK(a2 == b2); CHECK(aa == ba); CHECK(an == bn); CHECK(a4 == b4);
37+}
38+TEST_CASE("Metal rolling sum of squares and regression statistics match CPU", "[metal]") {
39+ auto m = metal(); if (!m) { SKIP("no Metal device"); } auto cpu = backend::make_cpu_backend();
40+ std::vector<f64> v; std::vector<i64> off; make_batch(3, 2000, v, off, 4); backend::Batch b{v, off};
41+ auto g = m->rolling_sum_squares(b, 22), c = cpu->rolling_sum_squares(b, 22);
42+ for (usize i = 0; i < v.size(); ++i) { if (std::isnan(c[i])) { CHECK(std::isnan(g[i])); continue; } CHECK_THAT(g[i], WithinRel(c[i], 1e-5)); }
43+ // regression: 5 series × 1000 rows × k=4
44+ Rng r(3); const usize S = 5, n = 1000, k = 4; std::vector<f64> X(S * n * k), y(S * n); std::vector<i64> ro(S + 1); for (usize s = 0; s <= S; ++s) ro[s] = static_cast<i64>(s * n);
45+ for (usize i = 0; i < S * n; ++i) { X[i * k] = 1.0; for (usize j = 1; j < k; ++j) X[i * k + j] = r.normal(); y[i] = 0.5 + X[i * k + 1] - 0.3 * X[i * k + 2] + r.normal(); }
46+ std::vector<f64> gx, gy, gyy, cx, cy, cyy; backend::KernelStats st; m->regression_stats(X, y, ro, k, gx, gy, gyy, &st); cpu->regression_stats(X, y, ro, k, cx, cy, cyy);
47+ for (usize i = 0; i < gx.size(); ++i) CHECK_THAT(gx[i], WithinAbs(cx[i], 1e-6 * std::max(1.0, std::abs(cx[i])) + 1e-9)); for (usize i = 0; i < gy.size(); ++i) CHECK_THAT(gy[i], WithinAbs(cy[i], 1e-6 * std::max(1.0, std::abs(cy[i])) + 1e-9)); for (usize s = 0; s < S; ++s) CHECK_THAT(gyy[s], WithinRel(cyy[s], 1e-6));
48+ // solve from GPU sufficient statistics and compare coefficients with the CPU QR solution
49+ Eigen::MatrixXd Xm(n, k); Eigen::VectorXd ym(n); for (usize i = 0; i < n; ++i) { ym[i] = y[i]; for (usize j = 0; j < k; ++j) Xm(i, j) = X[i * k + j]; }
50+ est::SufficientStats ss; ss.xtx = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(gx.data(), k, k); ss.xty = Eigen::Map<Eigen::VectorXd>(gy.data(), k); ss.yty = gyy[0]; ss.n = n;
51+ auto from_gpu = est::ols_from_stats(ss); auto ref = est::ols(Xm, ym); for (usize j = 0; j < k; ++j) CHECK_THAT(from_gpu.estimate.coefficients[j], WithinAbs(ref.estimate.coefficients[j], 1e-5));
52+}
53+TEST_CASE("Metal batched GARCH likelihood matches FP64 within tolerance", "[metal]") {
54+ auto m = metal(); if (!m) { SKIP("no Metal device"); } auto cpu = backend::make_cpu_backend();
55+ Rng rng(11); const usize S = 4, n = 2000; std::vector<f64> v(S * n); std::vector<i64> off(S + 1); for (usize s = 0; s <= S; ++s) off[s] = static_cast<i64>(s * n);
56+ for (usize s = 0; s < S; ++s) { f64 s2 = 0.04 / 0.05, e = 0; for (usize t = 0; t < n; ++t) { s2 = 0.04 + 0.05 * e * e + 0.9 * s2; e = std::sqrt(s2) * rng.normal(); v[s * n + t] = 0.01 + e; } }
57+ backend::Batch b{v, off};
58+ const usize P = 16; std::vector<f64> params(S * P * 4); for (usize s = 0; s < S; ++s) for (usize q = 0; q < P; ++q) { f64* th = params.data() + (s * P + q) * 4; th[0] = 0.01; th[1] = 0.02 + 0.004 * q; th[2] = 0.03 + 0.005 * q; th[3] = 0.85 + 0.005 * q; }
59+ backend::KernelStats st; auto g = m->garch_nll_batch(b, params, P, 4, 0, 0.0, &st); auto c = cpu->garch_nll_batch(b, params, P, 4, 0, 0.0);
60+ for (usize j = 0; j < S * P; ++j) { INFO("job " << j); CHECK_THAT(g[j], WithinRel(c[j], 2e-5)); }
61+ // argmin over the grid must agree
62+ for (usize s = 0; s < S; ++s) { usize ag = 0, ac = 0; for (usize q = 1; q < P; ++q) { if (g[s * P + q] < g[s * P + ag]) ag = q; if (c[s * P + q] < c[s * P + ac]) ac = q; } CHECK(ag == ac); }
63+ // GJR with gamma and Student-t fixed nu
64+ std::vector<f64> p5(S * P * 5); for (usize s = 0; s < S; ++s) for (usize q = 0; q < P; ++q) { f64* th = p5.data() + (s * P + q) * 5; th[0] = 0.01; th[1] = 0.03; th[2] = 0.03; th[3] = 0.02 + 0.002 * q; th[4] = 0.88; }
65+ auto g5 = m->garch_nll_batch(b, p5, P, 5, 1, 7.0, &st); auto c5 = cpu->garch_nll_batch(b, p5, P, 5, 1, 7.0); for (usize j = 0; j < S * P; ++j) CHECK_THAT(g5[j], WithinRel(c5[j], 5e-5));
66+}
67+TEST_CASE("Metal invalid parameters give +inf like the CPU", "[metal]") {
68+ auto m = metal(); if (!m) { SKIP("no Metal device"); }
69+ std::vector<f64> v(300, 0.01); std::vector<i64> off = {0, 300}; backend::Batch b{v, off}; std::vector<f64> bad = {0.0, 0.1, 0.6, 0.5};
70+ auto g = m->garch_nll_batch(b, bad, 1, 4, 0, 0.0); CHECK(std::isinf(g[0]));
71+}
72+#endif
73+
74+TEST_CASE("Backend resolution and fallback", "[backend]") {
75+ ExecutionPlan p; p.backend = Backend::cpu; std::string why; auto be = backend::resolve(p, &why); CHECK(be->kind() == Backend::cpu);
76+ p.backend = Backend::metal; p.allow_fallback = true; auto be2 = backend::resolve(p, &why); CHECK((be2->kind() == Backend::metal || be2->kind() == Backend::cpu));
77+ auto devs = backend::list_devices(); CHECK(!devs.empty()); CHECK(devs[0].kind == Backend::cpu); CHECK(devs[0].supports_fp64);
78+ for (auto& d : devs) if (d.kind == Backend::metal) CHECK_FALSE(d.supports_fp64);
79+}
80
Menu

Theme
↑↓ navigate ↵ open esc close / or ⌘K anywhere