// Metal backend host code (Objective-C++). See shaders/raid_kernels.metal for the kernels and the precision design.
#import <Metal/Metal.h>
#import <Foundation/Foundation.h>
#include "raid/backend/backend.hpp"
#include "raid/core/error.hpp"
#include "raid/core/log.hpp"
#include "raid/statistics/reductions.hpp"
#include <chrono>
#include <cmath>
#include <map>
#include <mutex>
#include <unistd.h>
#include <cstdlib>
#include <mach-o/dyld.h>

namespace raid { namespace metal { const char* msl_bundle(); } }

namespace raid::backend {

namespace {
struct Chunk { uint32_t seg, begin, end, pad; };
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(); } };
}

class MetalBackend final : public Backend {
 public:
  explicit MetalBackend(PrecisionPolicy policy) : policy_(policy) {
    @autoreleasepool {
      device_ = MTLCreateSystemDefaultDevice();
      if (!device_) { log::warn("metal: no device"); return; }
      queue_ = [device_ newCommandQueue];
      available_ = device_ != nil && queue_ != nil;
    }
  }
  ~MetalBackend() override = default;
  std::string name() const override { return "metal"; }
  raid::Backend kind() const override { return raid::Backend::metal; }
  bool available() const override { return available_; }
  Device device() const override {
    Device d; d.kind = raid::Backend::metal; d.supports_fp64 = false; d.unified_memory = device_.hasUnifiedMemory;
    d.name = device_ ? std::string([device_.name UTF8String]) : "none"; d.memory_bytes = device_.recommendedMaxWorkingSetSize; d.max_buffer_bytes = device_.maxBufferLength;
    d.compute_units = 0;  // Metal does not expose the GPU core count; recorded from system_profiler in the manifest
    d.details = std::string("families:") + ([device_ supportsFamily:MTLGPUFamilyApple9] ? "Apple9" : ([device_ supportsFamily:MTLGPUFamilyApple8] ? "Apple8" : "Apple7-")) + " metal3:" + ([device_ supportsFamily:MTLGPUFamilyMetal3] ? "yes" : "no") + " library:" + library_source_;
    return d;
  }

  // ---- library / pipelines ----------------------------------------------------------------------------------------
  bool ensure_library(KernelStats* st) {
    std::lock_guard<std::mutex> lk(m_);
    if (library_) return true;
    @autoreleasepool {
      NSError* err = nil; Timer t;
      // 1. prebuilt metallib: RAID_METALLIB env, compile-time path, next to the executable
      std::vector<std::string> candidates;
      const bool runtime_only = std::getenv("RAID_METAL_RUNTIME_ONLY") != nullptr;   // force runtime MSL compilation (benchmarks)
      if (const char* e = std::getenv("RAID_METALLIB")) candidates.push_back(e);
#ifdef RAID_METALLIB_PATH
      candidates.push_back(RAID_METALLIB_PATH);
#endif
      { 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"); } }
      for (auto& c : candidates) {
        if (runtime_only) break;
        if (access(c.c_str(), R_OK) != 0) continue;
        library_ = [device_ newLibraryWithURL:[NSURL fileURLWithPath:[NSString stringWithUTF8String:c.c_str()]] error:&err];
        if (library_) { library_source_ = "metallib:" + c; compile_seconds_ = t.s(); break; }
      }
      // 2. runtime compilation of the embedded MSL
      if (!library_) {
        MTLCompileOptions* opts = [MTLCompileOptions new]; opts.fastMathEnabled = NO;   // keep IEEE semantics (no reassociation): required by the double-float transforms
        // RAID_MSL_SALT appends a unique comment so that the system shader cache cannot serve a previous compilation
        // (used to measure genuinely cold runtime compilation in E7).
        std::string src = raid::metal::msl_bundle(); if (const char* salt = std::getenv("RAID_MSL_SALT")) src += "\n// salt " + std::string(salt) + "\n";
        library_ = [device_ newLibraryWithSource:[NSString stringWithUTF8String:src.c_str()] options:opts error:&err];
        if (!library_) { log::error("metal: shader compilation failed", {{"error", err ? std::string([[err localizedDescription] UTF8String]) : "?"}}); available_ = false; return false; }
        library_source_ = "runtime"; compile_seconds_ = t.s(); if (st) st->compiled_now = true;
      }
    }
    return true;
  }
  id<MTLComputePipelineState> pipeline(const std::string& fn) {
    std::lock_guard<std::mutex> lk(m_);
    auto it = pipelines_.find(fn); if (it != pipelines_.end()) return it->second;
    @autoreleasepool {
      NSError* err = nil; id<MTLFunction> f = [library_ newFunctionWithName:[NSString stringWithUTF8String:fn.c_str()]];
      RAID_REQUIRE(f != nil, unsupported, "metal: kernel not found: " + fn);
      id<MTLComputePipelineState> p = [device_ newComputePipelineStateWithFunction:f error:&err];
      RAID_REQUIRE(p != nil, unsupported, "metal: pipeline failed for " + fn + ": " + (err ? std::string([[err localizedDescription] UTF8String]) : "?"));
      pipelines_[fn] = p; return p;
    }
  }
  // ---- buffer pool (shared storage; unified memory) ---------------------------------------------------------------
  id<MTLBuffer> buffer(usize bytes, usize* allocs) {
    std::lock_guard<std::mutex> lk(m_);
    bytes = std::max<usize>(bytes, 16); usize cls = 1; while (cls < bytes) cls <<= 1;
    auto& v = pool_[cls]; if (!v.empty()) { auto b = v.back(); v.pop_back(); return b; }
    if (allocs) ++*allocs; ++total_allocs_;
    id<MTLBuffer> b = [device_ newBufferWithLength:cls options:MTLResourceStorageModeShared];
    RAID_REQUIRE(b != nil, memory, "metal: buffer allocation failed (" + std::to_string(cls) + " bytes)");
    return b;
  }
  void release(id<MTLBuffer> b) { std::lock_guard<std::mutex> lk(m_); pool_[b.length].push_back(b); }
  usize total_allocations() const { return total_allocs_; }

  // ---- helpers -----------------------------------------------------------------------------------------------------
  std::vector<Chunk> make_chunks(std::span<const i64> offsets, usize chunk_max, usize* n_segments) const {
    std::vector<Chunk> ch; const usize S = offsets.size() - 1; if (n_segments) *n_segments = S;
    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; }
      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}); }
    return ch;
  }
  usize chunk_size() const { return policy_.fp32_chunk_max == 0 ? 4096 : policy_.fp32_chunk_max; }
  static void finish(KernelStats* st, const char* k, const Timer& t, id<MTLCommandBuffer> cb, usize bin, usize bout, usize allocs, const std::string& note) {
    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;
  }
  // Fused realized kernel; returns 5 vectors combined in FP64 on the CPU (adds the chunk-boundary |x_i||x_{i-1}| terms).
  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) {
    Timer t; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
    usize S; auto chunks = make_chunks(b.offsets, chunk_size(), &S); const usize n = b.values.size(), C = chunks.size();
    s2.assign(S, 0); adj.assign(S, 0); neg.assign(S, 0); pos.assign(S, 0); s4.assign(S, 0);
    if (n == 0 || C == 0) { finish(st, "realized_fused", t, nil, 0, 0, 0, library_source_); return; }
    @autoreleasepool {
      id<MTLBuffer> bx = buffer(n * 4, &allocs), bc = buffer(C * sizeof(Chunk), &allocs), bo = buffer(C * 5 * 8, &allocs);
      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();
      std::memcpy(bc.contents, chunks.data(), C * sizeof(Chunk));
      id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder];
      const uint32_t op = single_op < 0 ? 0 : static_cast<uint32_t>(single_op);
      id<MTLComputePipelineState> p = pipeline(single_op < 0 ? "realized_fused" : "segment_reduce");
      [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];
      [enc dispatchThreadgroups:MTLSizeMake(C, 1, 1) threadsPerThreadgroup:MTLSizeMake(256, 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
      RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
      const float* o = static_cast<const float*>(bo.contents);
      for (usize c = 0; c < C; ++c) { const usize s = chunks[c].seg;
        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]; }
        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; }
        // chunk boundary term for the adjacent product (FP64 on the CPU)
        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]); }
      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));
      release(bx); release(bc); release(bo);
    }
  }

  // ---- Backend interface ---------------------------------------------------------------------------------------------
  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; }
  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; }
  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); }
  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; }
  /// All five realized accumulators in a single fused pass (used by est::realized_measures_batch when the backend is Metal).
  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); }

  std::vector<f64> rolling_sum_squares(const Batch& b, usize window, KernelStats* st) override {
    Timer t; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
    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;
    RAID_REQUIRE(window <= 65536, unsupported, "metal rolling window too large for the direct kernel");
    @autoreleasepool {
      id<MTLBuffer> bx = buffer(n * 4, &allocs), bb = buffer(S * 4, &allocs), bs = buffer(n * 4, &allocs), bo = buffer(n * 4, &allocs);
      float* xf = static_cast<float*>(bx.contents); for (usize i = 0; i < n; ++i) xf[i] = static_cast<float>(b.values[i]);
      uint32_t* sb = static_cast<uint32_t*>(bb.contents); uint32_t* so = static_cast<uint32_t*>(bs.contents);
      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); }
      const uint32_t w = static_cast<uint32_t>(window), nn = static_cast<uint32_t>(n);
      id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder]; auto p = pipeline("rolling_sumsq");
      [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];
      [enc dispatchThreads:MTLSizeMake(n, 1, 1) threadsPerThreadgroup:MTLSizeMake(std::min<usize>(256, p.maxTotalThreadsPerThreadgroup), 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
      RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
      const float* o = static_cast<const float*>(bo.contents); for (usize i = 0; i < n; ++i) out[i] = o[i];
      finish(st, "rolling_sumsq", t, cb, n * 8, n * 8, allocs, library_source_); release(bx); release(bb); release(bs); release(bo);
    }
    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; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
    RAID_REQUIRE(k >= 1 && k <= 8, unsupported, "metal regression_stats supports k <= 8");
    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;
    xtx.assign(S * k * k, 0); xty.assign(S * k, 0); yty.assign(S, 0); if (rows == 0) return;
    @autoreleasepool {
      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);
      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]);
      std::memcpy(bc.contents, chunks.data(), C * sizeof(Chunk)); const uint32_t kk = static_cast<uint32_t>(k);
      id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder]; auto p = pipeline("regression_stats");
      [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];
      [enc dispatchThreadgroups:MTLSizeMake(C, 1, 1) threadsPerThreadgroup:MTLSizeMake(256, 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
      RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
      const float* o = static_cast<const float*>(bo.contents);
      for (usize c = 0; c < C; ++c) { const usize s = chunks[c].seg; usize q = 0;
        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; }
        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; }
        yty[s] += static_cast<f64>(o[(c * nstat + q) * 2]) + o[(c * nstat + q) * 2 + 1]; }
      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);
    }
  }

  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 {
    Timer t; usize allocs = 0; RAID_REQUIRE(available_ && ensure_library(st), unsupported, "metal unavailable");
    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;
    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)");
    @autoreleasepool {
      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);
      float* rf = static_cast<float*>(br.contents);
      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);
      float* pf = static_cast<float*>(bp.contents); uint32_t* jb = static_cast<uint32_t*>(bj.contents);
      const bool gjr = model_code == 1, eg = model_code == 2;
      // Per-series scaling to unit variance (FP32 conditioning): r' = r/sd, μ' = μ/sd, ω' = ω/sd² (GARCH/GJR) or ω' = ω + (1−β)·log(1/sd²) (EGARCH);
      // the NLL of the original data is NLL' + n·log(sd). Scales are chosen in FP64 on the host.
      std::vector<f64> sd(S, 1.0);
      for (usize s = 0; s < S; ++s) { auto x = r.series(s); f64 m = 0; for (f64 v : x) m += v; m /= std::max<usize>(1, x.size()); f64 var = 0; for (f64 v : x) var += (v - m) * (v - m); var /= std::max<usize>(1, x.size()); sd[s] = var > 0 ? std::sqrt(var) : 1.0; for (i64 i = r.offsets[s]; i < r.offsets[s + 1]; ++i) rf[i] = static_cast<float>(r.values[static_cast<usize>(i)] / sd[s]); }
      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]);
        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;
          // 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
          const f64 beta_s = th[(gjr || eg) ? 4 : 3];
          pf[j * 8 + 0] = static_cast<float>(th[0] / sd[s]); pf[j * 8 + 1] = static_cast<float>(eg ? th[1] + (1.0 - beta_s) * std::log(1.0 / (sd[s] * sd[s])) : th[1] / (sd[s] * sd[s])); pf[j * 8 + 2] = static_cast<float>(th[2]);
          pf[j * 8 + 3] = (gjr || eg) ? static_cast<float>(th[3]) : 0.0f; pf[j * 8 + 4] = static_cast<float>(beta_s);
          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);
          // FP64 density constants (Metal has no lgamma): kconst and E|z|
          const f64 pi = 3.14159265358979323846;
          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))); }
          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)); }
          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; }
        // initial variance: sample variance of (r − mu) with mu from the first parameter vector of the series (FP64)
        const f64 mu = params[(s * n_params) * p] / sd[s]; auto x = r.series(s); f64 m1 = 0; for (f64 v : x) m1 += v / sd[s] - mu; m1 /= static_cast<f64>(x.size()); f64 var0 = 0; for (f64 v : x) var0 += (v / sd[s] - mu - m1) * (v / sd[s] - mu - m1); var0 /= static_cast<f64>(x.size()); v0[s] = static_cast<float>(var0); }
      const uint32_t model = static_cast<uint32_t>(model_code), nj = static_cast<uint32_t>(J);
      id<MTLCommandBuffer> cb = [queue_ commandBuffer]; id<MTLComputeCommandEncoder> enc = [cb computeCommandEncoder]; auto pl = pipeline("garch_nll");
      [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];
      [enc dispatchThreads:MTLSizeMake(J, 1, 1) threadsPerThreadgroup:MTLSizeMake(std::min<usize>(64, pl.maxTotalThreadsPerThreadgroup), 1, 1)]; [enc endEncoding]; [cb commit]; [cb waitUntilCompleted];
      RAID_REQUIRE(cb.status == MTLCommandBufferStatusCompleted, numerical, "metal: command buffer failed");
      const float* o = static_cast<const float*>(bo.contents); for (usize j = 0; j < J; ++j) { const usize s = j / n_params; out[j] = std::isfinite(o[j]) ? static_cast<f64>(o[j]) + static_cast<f64>(r.size(s)) * std::log(sd[s]) : o[j]; }
      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);
    }
    return out;
  }
  f64 compile_seconds() const { return compile_seconds_; }
  std::string library_source() const { return library_source_; }

 private:
  PrecisionPolicy policy_; id<MTLDevice> device_ = nil; id<MTLCommandQueue> queue_ = nil; id<MTLLibrary> library_ = nil; bool available_{false};
  std::map<std::string, id<MTLComputePipelineState>> pipelines_; std::map<usize, std::vector<id<MTLBuffer>>> pool_; std::mutex m_; usize total_allocs_{0};
  std::string library_source_{"none"}; f64 compile_seconds_{0};
};

std::shared_ptr<Backend> make_metal_backend(const PrecisionPolicy& policy) {
  static std::mutex m; static std::shared_ptr<MetalBackend> cached; std::lock_guard<std::mutex> lk(m);
  if (cached) return cached->available() ? cached : nullptr;   // one device/queue/library/pipeline cache/buffer pool per process
  cached = std::make_shared<MetalBackend>(policy); return cached->available() ? cached : nullptr;
}

/// Fused realized accumulators when the backend is Metal (falls back to the four generic calls otherwise).
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) {
  if (auto* m = dynamic_cast<MetalBackend*>(&be)) { m->realized_fused(b, s2, adj, neg, pos, s4, st); return true; }
  return false;
}
f64 metal_compile_seconds(Backend& be) { if (auto* m = dynamic_cast<MetalBackend*>(&be)) return m->compile_seconds(); return 0; }

}  // namespace raid::backend