Merge pull request 'engram: make the ggml batch-cosine strategy actually compute in fp32 (recall 0.9933 -> 0.9987)' (#121) from improve/ggml-cosine-fp32-and-init into feat/engram-ggml-cosine-batch

#116 shipped the ggml strategy at 0.9933 id-recall against the CPU oracle
while the hand-rolled Metal kernel it replaced scored 0.9997 — a ~150x worse
error margin. That was not an inherent property of ggml. It was a usage bug in
this file, and this commit fixes it.

ggml-metal has two F32xF32 matmul kernels and picks between them purely on
ne11, the number of B rows, which for us is the query-batch size:

  ne11 <= 8  -> kernel_mul_mv_ext_f32_f32_* / kernel_mul_mv_f32_f32_*,
                templated <float, float> — genuine F32.
  ne11 >  8  -> kernel_mul_mm_f32_f32, templated
                <half, half4x4, simdgroup_half8x8, half, half2x4, ...> —
                BOTH operands narrowed to F16, despite F32 tensors on both
                sides.

The old code issued one ggml_mul_mat with ne11 = nq (300 in the benchmark),
landing squarely on the F16 path. The file's own header comment asserted the
opposite ("computes in F32 on the Metal backend"); that claim was wrong and is
replaced with the measurement.

Fix: emit ceil(nq/8) mul_mats over ne11<=8 ggml_view_2d slices of one query
tensor, all expanded into ONE graph and one ggml_backend_graph_compute, so the
node matrix is still uploaded and shared exactly once. EL_GGML_MULMAT_CHUNK
overrides the 8; setting it >= nq reproduces the old behaviour exactly, which
is also how the before/after below was measured in a single binary.

Measured, real store snapshot, 13415 live embedded nodes, dim=768, 300 real
queries, vs the CPU double-accumulated oracle (vindex_bench, offline copy of
the store — no live service touched):

  id-recall   same-rank |Δdist| max   mean
  old (ne11=300)   0.9933   6.80e-05   1.43e-05
  new (ne11<=8)    0.9987   4.77e-07   9.30e-08
  hand-rolled      0.9997   3.58e-07   7.55e-08

~145x better max error, ~154x better mean — now the same order of magnitude as
the hand-rolled kernel rather than 150x off it.

The cost is real and is documented rather than buried. Median of 15 reps of
the whole batch_multi() call, three runs: 13.2-14.4ms unchunked, 19.9-20.2ms
chunked, 17.7-18.0ms hand-rolled. Correctness costs ~+6.7ms per 300-query
batch and leaves ggml ~12% behind the hand-rolled kernel instead of ~35%
ahead. It cannot be recovered inside ggml: an fp32 matmul on Metal must
re-stream the node matrix once per <=8 queries, and ggml's Metal backend ships
no fp32 TILED matmul, so "fast" and "fp32" are genuinely exclusive there.

Two things that did NOT work, recorded so nobody retries them:

  - ggml_mul_mat_set_prec(t, GGML_PREC_F32) does nothing here. Error was
    bit-identical with and without it (1.038e-05 either way) — ggml-metal has
    no F32-accumulating mul_mm kernel to switch to. ne11 is the only lever.
  - The ACCEL/BLAS device looked excellent in an isolated compute-only probe
    (3.4-4.0ms, mean |Δdot| 1.5e-08) but is dominated on BOTH axes end-to-end
    (0.191 ms/query at 0.9973 recall vs 0.125-0.142 at 0.9987), because the
    probe was not competing for the same CPU cores the real call path is. It
    stays reachable via EL_GGML_DEVICE as a no-Metal fallback, labelled as
    measured-and-rejected, not as a recommendation.

Also corrected: the ~7.8s "cold start" blamed on this file is not this file
re-initialising per call — init was already cached. It is Apple's shader cache
missing on ggml's embedded metallib (~650 kernels), keyed on the library and
shared across processes: the first load on a machine reports
"loaded in 7.670 sec", the next run of a *different* binary reports 0.009 sec.
Once per machine per ggml version, not once per process, and not ours to fix.
Warm ggml init is 44-53ms vs 36-117ms for the hand-rolled strategy.

Loading only libggml-metal.so instead of every plugin in the directory is kept
for tidiness, and explicitly documented as NOT a speedup: 44.7-52.4ms against
46.9-58.9ms, the same number inside noise.

The -2.0 sentinel contract is unchanged and re-verified at batch sizes that
straddle the chunk boundary (1,7,8,9,16,17,33), plus NULL rows, dim
mismatches, zero-norm rows, and an all-invalid population. Notably the old
ne11=300 path fails that same check at a 2e-6 cosine tolerance with 2299
mismatches, which is an independent confirmation of the defect.

lang/runtime/build_vindex_bench.sh

@@ -0,0 +1,51 @@
+#!/bin/bash
+# build_vindex_bench.sh — build the vindex_bench oracle/proof harness, with
+# the real ggml + hand-rolled-Metal batch-cosine strategies on Darwin and a
+# zero-dependency CPU-only stub everywhere else. Mirrors the two-step recipe
+# documented in vindex_bench.c's own header comment; this script exists so
+# that recipe is one command, not a copy-pasted paragraph.
+#
+# Darwin build links FOUR strategy translation units:
+#   eg_cosine_batch.c                    — the Factory (always)
+#   eg_cosine_batch_strategy_cpu.c       — universal fallback (always)
+#   eg_cosine_batch_strategy_ggml.c      — ggml + dynamic Metal backend plugin
+#   eg_cosine_batch_strategy_metal_hand.m — PR #114's original hand-rolled
+#                                            Metal shader, preserved as one
+#                                            selectable strategy
+# plus -DEG_HAVE_STRATEGY_GGML -DEG_HAVE_STRATEGY_METAL_HAND so the Factory
+# (and vindex_bench.c's own direct strategy comparison) knows both exist.
+#
+# ggml is resolved via `brew --prefix ggml` when available (portable across
+# Intel /usr/local and Apple Silicon /opt/homebrew installs), falling back to
+# /opt/homebrew if brew isn't on PATH. Override with GGML_PREFIX=... env var.
+#
+# Usage: ./build_vindex_bench.sh [output_path]
+set -euo pipefail
+cd "$(dirname "$0")"
+OUT="${1:-./vindex_bench}"
+CC="${CC:-cc}"
+
+if [ "$(uname -s)" = "Darwin" ]; then
+  GGML_PREFIX="${GGML_PREFIX:-$(brew --prefix ggml 2>/dev/null || echo /opt/homebrew)}"
+  echo "== Darwin: building with the ggml + hand-rolled-Metal strategies (ggml prefix: $GGML_PREFIX) =="
+
+  "$CC" -O2 -std=c11 -x objective-c \
+    -c eg_cosine_batch_strategy_metal_hand.m -o /tmp/eg_cosine_batch_strategy_metal_hand.o \
+    -framework Metal -framework Foundation
+
+  "$CC" -O2 -std=c11 -DEG_HAVE_STRATEGY_GGML -DEG_HAVE_STRATEGY_METAL_HAND -w \
+    -I"$GGML_PREFIX/include" \
+    vindex_bench.c engram_vindex.c \
+    eg_cosine_batch.c eg_cosine_batch_strategy_cpu.c eg_cosine_batch_strategy_ggml.c \
+    /tmp/eg_cosine_batch_strategy_metal_hand.o \
+    -L"$GGML_PREFIX/lib" -lggml -lggml-base \
+    -Wl,-rpath,"$GGML_PREFIX/lib" \
+    -lm -framework Metal -framework Foundation -o "$OUT"
+else
+  echo "== non-Darwin: building with the CPU-only fallback strategy (no ggml, no Metal) =="
+  "$CC" -O2 -std=c11 -w vindex_bench.c engram_vindex.c \
+    eg_cosine_batch.c eg_cosine_batch_strategy_cpu.c \
+    -lm -o "$OUT"
+fi
+
+echo "built: $OUT"

lang/runtime/eg_cosine_batch.c

@@ -0,0 +1,121 @@
+/* eg_cosine_batch.c — the Factory. Implements the stable public interface
+ * declared in eg_cosine_batch.h by selecting ONE concrete
+ * EgCosineBatchStrategy (eg_cosine_batch_strategy.h) and dispatching every
+ * call to it. This is the ONLY file that branches on EG_HAVE_STRATEGY_*
+ * (build-time: which strategy .c/.m files were actually compiled in for
+ * this platform) — call sites never see those macros.
+ *
+ * Selection is lazy (first call) and cached — mirrors the lazy-init caching
+ * every individual strategy already does internally, so there is no added
+ * per-call cost after the first.
+ *
+ * Selection mechanism (env var + build-time + runtime capability probe, all
+ * three, exactly as directed):
+ *   - BUILD-TIME decides which strategies exist to choose from at all: a
+ *     Darwin build compiles+links the ggml strategy and the hand-rolled
+ *     Metal strategy (EG_HAVE_STRATEGY_GGML / EG_HAVE_STRATEGY_METAL_HAND
+ *     both defined); a non-Darwin build compiles neither, matching PR #114's
+ *     original Linux behavior exactly (CPU-fallback only, no Objective-C
+ *     compiler or Metal frameworks required).
+ *   - RUNTIME CAPABILITY PROBE: each candidate strategy's own available()
+ *     does the real, cheap-after-first-call check (device present, backend
+ *     plugin loaded, pipeline compiled) — never assumed from build-time
+ *     alone. A build that HAS the ggml strategy compiled in but is running
+ *     on hardware/software where it can't actually initialize (backend
+ *     plugin missing, no GPU) correctly falls through to the next candidate.
+ *   - ENV VAR gives explicit, debuggable override for either axis:
+ *       EL_COSINE_BATCH_STRATEGY = "ggml" | "metal" | "cpu" | unset/"auto"
+ *         forces a specific strategy (falling back to cpu if the forced one
+ *         isn't actually available), or leaves the default auto-preference
+ *         order in place.
+ *       EL_METAL_COSINE = 0/n/N/f/F   (back-compat with PR #114's vindex_bench
+ *         gate) disables ALL GPU-backed strategies outright, same as before.
+ *
+ * DEFAULT preference order when nothing is forced: ggml, then hand-rolled
+ * Metal, then CPU fallback — first candidate whose available() reports true
+ * wins. This is what makes "stop hand-rolling GPU kernels, use ggml" real
+ * rather than nominal: ggml is what actually runs by default on this
+ * machine today (see the PR body for the measured numbers backing that).
+ */
+#include "eg_cosine_batch.h"
+#include "eg_cosine_batch_strategy.h"
+
+#include <stdlib.h>
+#include <string.h>
+
+static bool g_selected = false;
+static const EgCosineBatchStrategy* g_active = NULL;
+
+static bool eg_env_truthy_off(const char* v) {
+    return v && (v[0]=='0' || v[0]=='n' || v[0]=='N' || v[0]=='f' || v[0]=='F');
+}
+
+static const EgCosineBatchStrategy* eg_select_strategy(void) {
+    if (g_selected) return g_active;
+    g_selected = true;
+
+    const EgCosineBatchStrategy* cpu = eg_cosine_batch_strategy_cpu();
+    const char* force = getenv("EL_COSINE_BATCH_STRATEGY");
+    const char* legacy_off = getenv("EL_METAL_COSINE");
+
+    if (eg_env_truthy_off(legacy_off)) { g_active = cpu; return g_active; }
+
+    if (force && strcmp(force, "cpu") == 0) { g_active = cpu; return g_active; }
+
+    if (force && strcmp(force, "ggml") == 0) {
+#ifdef EG_HAVE_STRATEGY_GGML
+        const EgCosineBatchStrategy* s = eg_cosine_batch_strategy_ggml();
+        if (s->available()) { g_active = s; return g_active; }
+#endif
+        g_active = cpu; return g_active;
+    }
+
+    if (force && strcmp(force, "metal") == 0) {
+#ifdef EG_HAVE_STRATEGY_METAL_HAND
+        const EgCosineBatchStrategy* s = eg_cosine_batch_strategy_metal_hand();
+        if (s->available()) { g_active = s; return g_active; }
+#endif
+        g_active = cpu; return g_active;
+    }
+
+    /* auto (unset, or any other value): ggml -> metal-hand -> cpu, first
+     * available wins. */
+#ifdef EG_HAVE_STRATEGY_GGML
+    {
+        const EgCosineBatchStrategy* s = eg_cosine_batch_strategy_ggml();
+        if (s->available()) { g_active = s; return g_active; }
+    }
+#endif
+#ifdef EG_HAVE_STRATEGY_METAL_HAND
+    {
+        const EgCosineBatchStrategy* s = eg_cosine_batch_strategy_metal_hand();
+        if (s->available()) { g_active = s; return g_active; }
+    }
+#endif
+    g_active = cpu;
+    return g_active;
+}
+
+bool eg_cosine_batch_available(void) {
+    return eg_select_strategy()->available();
+}
+
+const char* eg_cosine_batch_strategy_name(void) {
+    return eg_select_strategy()->name;
+}
+
+bool eg_cosine_batch(const float* query, int32_t qdim,
+                      const float* const* node_ptrs,
+                      const int32_t* node_dims,
+                      int32_t n,
+                      double* out_scores) {
+    return eg_select_strategy()->batch(query, qdim, node_ptrs, node_dims, n, out_scores);
+}
+
+bool eg_cosine_batch_multi(const float* queries, int32_t qdim, int32_t nq,
+                            const float* const* node_ptrs,
+                            const int32_t* node_dims,
+                            int32_t n,
+                            double* out_scores) {
+    return eg_select_strategy()->batch_multi(queries, qdim, nq, node_ptrs, node_dims, n, out_scores);
+}

lang/runtime/eg_cosine_batch.h

@@ -0,0 +1,106 @@
+/* eg_cosine_batch.h — stable Adapter interface over batch-cosine-similarity
+ * BACKEND STRATEGIES. This header is the ONE thing call sites (el_runtime.c,
+ * vindex_bench.c, ...) talk to. Plain C11, safe to #include on every
+ * platform — the symbols declared here always exist and always link,
+ * regardless of what backend actually runs underneath. Zero #ifdef at call
+ * sites: which concrete strategy executes (ggml/Metal, hand-rolled Metal, or
+ * the always-false CPU fallback) is resolved once, lazily, inside
+ * eg_cosine_batch.c's factory — see eg_cosine_batch_strategy.h for that.
+ *
+ * This supersedes eg_metal_cosine.h (PR #114's single hand-rolled-Metal-only
+ * bridge). The contract is UNCHANGED — same shapes, same sentinel, same
+ * never-partial guarantee, same "caller must always be prepared to fall back
+ * to its own scalar per-node loop" rule — only the name changed, because the
+ * thing behind it is no longer "the Metal bridge," it is "whichever batch-
+ * cosine strategy the factory picked." eg_metal_cosine.h's original doc
+ * comments (byte-for-byte, this file is the direct descendant) are preserved
+ * below since they remain the precise spec any strategy must honor.
+ *
+ * On ANY failure at ANY step — no compute device, compile/init error, alloc
+ * failure, bad args — every function here returns false and writes nothing.
+ * Out-params are either fully populated or left completely untouched, never
+ * partial. Callers MUST always be prepared to fall back to their own scalar
+ * per-node CPU loop unconditionally. These functions must never crash, throw,
+ * or hang the calling process — several call sites run inside a long-lived
+ * daemon's request-handling hot path.
+ */
+#ifndef EG_COSINE_BATCH_H
+#define EG_COSINE_BATCH_H
+
+#include <stdint.h>
+#include <stdbool.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Batched cosine similarity: one query vector against `n` node vectors.
+ *
+ *   query      — qdim floats, the query embedding. Raw/unnormalized.
+ *   qdim       — query dimensionality (e.g. 768 for nomic-embed-text).
+ *   node_ptrs  — array of n pointers, node_ptrs[i] pointing at a (possibly
+ *                differently-owned, possibly NULL) float vector for node i.
+ *                NOT required to be contiguous — every strategy performs the
+ *                gather into a packed row-major matrix internally, exactly
+ *                mirroring how EngramNode.emb is one malloc per node.
+ *   node_dims  — array of n ints, node_dims[i] = that node's real emb_dim
+ *                (0 or mismatched vs qdim => that node scores -2.0, matching
+ *                eg_cosine's null/dim-mismatch/zero-norm sentinel exactly).
+ *   n          — number of nodes.
+ *   out_scores — caller-owned array of n doubles; out_scores[i] is filled
+ *                with the cosine similarity of node i against query, or
+ *                -2.0 for a null/dim-mismatched/zero-norm node — bit-for-bit
+ *                the same contract as eg_cosine(node_ptrs[i], query, qdim).
+ *
+ * Returns true iff a real backend strategy ran and out_scores was fully
+ * populated. Returns false (out_scores left untouched) on ANY failure or
+ * unavailability — no compute device, compile/init failure, allocation
+ * failure, n<=0, qdim<=0, null query/node_ptrs/node_dims/out_scores.
+ */
+bool eg_cosine_batch(const float* query, int32_t qdim,
+                      const float* const* node_ptrs,
+                      const int32_t* node_dims,
+                      int32_t n,
+                      double* out_scores);
+
+/* True iff a real (non-CPU-fallback) strategy is available right now (cheap
+ * after the first call — cached). Purely informational (e.g. a startup log
+ * line or /api/stats field); callers should still treat a false return from
+ * eg_cosine_batch()/eg_cosine_batch_multi() itself as the authoritative
+ * fallback signal, not this function. */
+bool eg_cosine_batch_available(void);
+
+/* Which concrete strategy is currently selected — "ggml", "metal-hand",
+ * or "cpu-fallback". Purely informational/diagnostic, same spirit as
+ * eg_cosine_batch_available(). Never NULL. */
+const char* eg_cosine_batch_strategy_name(void);
+
+/* Multi-query batched cosine: nq query vectors against the SAME n node
+ * vectors, in one call. A real strategy uploads/prepares the node population
+ * ONCE and reuses it for every query, instead of nq separate
+ * eg_cosine_batch() calls each paying the full gather+upload cost — PR #114
+ * measured this necessary: at N~=13.7k/dim=768, repeating the single-query
+ * call per query was slower than the CPU baseline; batching queries together
+ * is what makes a GPU-backed path a real win at this shape. Use this
+ * whenever multiple queries will run against an unchanged (or
+ * rarely-changing) node population; use eg_cosine_batch() for a genuinely
+ * one-off comparison.
+ *
+ *   queries    — nq*qdim floats, row-major (query i at queries+i*qdim).
+ *   out_scores — caller-owned nq*n doubles, row-major
+ *                (out_scores[i*n+j] = cosine(queries[i], node j)), same
+ *                -2.0 sentinel semantics as eg_cosine_batch().
+ *
+ * Returns true iff a real strategy ran and out_scores was fully populated
+ * (all nq*n entries); false (untouched) on any failure/unavailability. */
+bool eg_cosine_batch_multi(const float* queries, int32_t qdim, int32_t nq,
+                            const float* const* node_ptrs,
+                            const int32_t* node_dims,
+                            int32_t n,
+                            double* out_scores);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* EG_COSINE_BATCH_H */

lang/runtime/eg_cosine_batch.metal

@@ -0,0 +1,156 @@
+/* eg_cosine_batch.metal — batched cosine similarity, one query vs N node vectors.
+ *
+ * GPU-shaped counterpart to eg_cosine() in el_runtime.c: same math, same
+ * dim-mismatch/zero-norm sentinel (-2.0), applied to N independent rows in
+ * parallel instead of one pair at a time in a CPU loop.
+ *
+ * Semantics MUST match eg_cosine() exactly:
+ *   - inputs are raw, UNNORMALIZED vectors (nomic-embed-text magnitudes are
+ *     not 1.0) — this kernel computes the full dot/(|a|*|b|) cosine, not a
+ *     plain dot product.
+ *   - a node whose declared dim differs from the query dim, or whose norm is
+ *     zero, scores exactly -2.0 (below any valid cosine in [-1,1]), so a
+ *     caller doing `if (score > threshold)` behaves identically whether the
+ *     scalar or the batched path filled the array.
+ *
+ * Precision: Apple GPUs do not support double in Metal Shading Language —
+ * everything here is float32. eg_cosine accumulates in CPU double, but its
+ * *inputs* are float32 embeddings, so the achievable precision ceiling is
+ * bounded by the input data regardless of accumulator width. To keep the
+ * float32 reduction from drifting relative to the double-accumulated CPU
+ * result across dim=768 terms, each thread accumulates with 4 independent
+ * partial sums (unrolled) rather than one running scalar — the same
+ * error-reduction trick already used by the CPU brute-force loop in
+ * vindex_bench.c. The measured float-vs-double delta is reported in the PR
+ * description; this is not assumed to be "close enough" without measurement.
+ */
+#include <metal_stdlib>
+using namespace metal;
+
+/* Per-dispatch invariants. `dim` is the query's dimensionality — the
+ * dimensionality every comparable node vector must match. */
+struct EgCosineParams {
+    uint n;      /* number of node rows */
+    uint dim;    /* vector width (both query and node rows are `dim` wide in
+                  * the packed buffer; node_dims[] carries each node's REAL
+                  * embedded dim for the mismatch check) */
+};
+
+/* One thread per node row. node_matrix is n*dim floats, row-major, packed at
+ * `dim` stride regardless of a row's real dim (the CPU side zero-pads or
+ * skips packing rows that don't match — see eg_cosine_batch_metal in
+ * eg_metal_cosine.m for the exact packing contract). node_dims[i] is the
+ * node's true emb_dim, used only for the mismatch sentinel — never used to
+ * index, since every row is packed at uniform `dim` stride. */
+kernel void eg_cosine_batch_kernel(
+    device const float*  query        [[buffer(0)]],
+    device const float*  node_matrix  [[buffer(1)]],
+    device const int*    node_dims    [[buffer(2)]],
+    constant EgCosineParams& p        [[buffer(3)]],
+    device float*        out_scores   [[buffer(4)]],
+    uint gid [[thread_position_in_grid]])
+{
+    if (gid >= p.n) return;
+
+    if (node_dims[gid] != int(p.dim)) {
+        out_scores[gid] = -2.0f;
+        return;
+    }
+
+    device const float* row = node_matrix + (uint64_t)gid * (uint64_t)p.dim;
+
+    /* 4-way partial accumulation — same shape as vindex_bench.c's brute_topk
+     * unroll, done here for float32 accuracy rather than raw throughput. */
+    float dot0 = 0.0f, dot1 = 0.0f, dot2 = 0.0f, dot3 = 0.0f;
+    float na0  = 0.0f, na1  = 0.0f, na2  = 0.0f, na3  = 0.0f;
+    float nb0  = 0.0f, nb1  = 0.0f, nb2  = 0.0f, nb3  = 0.0f;
+
+    uint d = 0;
+    uint dim4 = p.dim & ~3u;
+    for (; d < dim4; d += 4) {
+        float a0 = row[d],   b0 = query[d];
+        float a1 = row[d+1], b1 = query[d+1];
+        float a2 = row[d+2], b2 = query[d+2];
+        float a3 = row[d+3], b3 = query[d+3];
+        dot0 += a0*b0; dot1 += a1*b1; dot2 += a2*b2; dot3 += a3*b3;
+        na0  += a0*a0; na1  += a1*a1; na2  += a2*a2; na3  += a3*a3;
+        nb0  += b0*b0; nb1  += b1*b1; nb2  += b2*b2; nb3  += b3*b3;
+    }
+    float dot = (dot0 + dot1) + (dot2 + dot3);
+    float na  = (na0  + na1)  + (na2  + na3);
+    float nb  = (nb0  + nb1)  + (nb2  + nb3);
+    for (; d < p.dim; d++) {
+        float a = row[d], b = query[d];
+        dot += a*b; na += a*a; nb += b*b;
+    }
+
+    if (na <= 0.0f || nb <= 0.0f) {
+        out_scores[gid] = -2.0f;
+        return;
+    }
+    out_scores[gid] = dot / sqrt(na * nb);
+}
+
+/* ── multi-query variant ──────────────────────────────────────────────────
+ * Same per-pair math as eg_cosine_batch_kernel, but amortizes ONE upload of
+ * node_matrix (the expensive part at real store size — 13k*768 floats is
+ * ~42MB) across `nq` queries instead of re-uploading it once per query.
+ * Measured need: a naive one-query-at-a-time loop calling the single-query
+ * kernel nq times was SLOWER than the CPU oracle at N≈13.7k (re-gather +
+ * re-upload dominated the actual compute) — this is the fix, not a
+ * hypothetical optimization.
+ *
+ * 2D grid: x = node index [0,n), y = query index [0,nq). out_scores is
+ * nq*n, row-major by query (out_scores[qid*n + nid]). */
+struct EgCosineMultiParams { uint n; uint dim; uint nq; };
+
+kernel void eg_cosine_batch_multi_kernel(
+    device const float*  queries      [[buffer(0)]],   /* nq*dim */
+    device const float*  node_matrix  [[buffer(1)]],   /* n*dim */
+    device const int*    node_dims    [[buffer(2)]],   /* n */
+    constant EgCosineMultiParams& p   [[buffer(3)]],
+    device float*        out_scores   [[buffer(4)]],   /* nq*n */
+    uint2 gid [[thread_position_in_grid]])
+{
+    uint nid = gid.x, qid = gid.y;
+    if (nid >= p.n || qid >= p.nq) return;
+
+    uint64_t out_idx = (uint64_t)qid * (uint64_t)p.n + (uint64_t)nid;
+
+    if (node_dims[nid] != int(p.dim)) {
+        out_scores[out_idx] = -2.0f;
+        return;
+    }
+
+    device const float* row   = node_matrix + (uint64_t)nid * (uint64_t)p.dim;
+    device const float* query = queries     + (uint64_t)qid * (uint64_t)p.dim;
+
+    float dot0 = 0.0f, dot1 = 0.0f, dot2 = 0.0f, dot3 = 0.0f;
+    float na0  = 0.0f, na1  = 0.0f, na2  = 0.0f, na3  = 0.0f;
+    float nb0  = 0.0f, nb1  = 0.0f, nb2  = 0.0f, nb3  = 0.0f;
+
+    uint d = 0;
+    uint dim4 = p.dim & ~3u;
+    for (; d < dim4; d += 4) {
+        float a0 = row[d],   b0 = query[d];
+        float a1 = row[d+1], b1 = query[d+1];
+        float a2 = row[d+2], b2 = query[d+2];
+        float a3 = row[d+3], b3 = query[d+3];
+        dot0 += a0*b0; dot1 += a1*b1; dot2 += a2*b2; dot3 += a3*b3;
+        na0  += a0*a0; na1  += a1*a1; na2  += a2*a2; na3  += a3*a3;
+        nb0  += b0*b0; nb1  += b1*b1; nb2  += b2*b2; nb3  += b3*b3;
+    }
+    float dot = (dot0 + dot1) + (dot2 + dot3);
+    float na  = (na0  + na1)  + (na2  + na3);
+    float nb  = (nb0  + nb1)  + (nb2  + nb3);
+    for (; d < p.dim; d++) {
+        float a = row[d], b = query[d];
+        dot += a*b; na += a*a; nb += b*b;
+    }
+
+    if (na <= 0.0f || nb <= 0.0f) {
+        out_scores[out_idx] = -2.0f;
+        return;
+    }
+    out_scores[out_idx] = dot / sqrt(na * nb);
+}

lang/runtime/eg_cosine_batch_strategy.h

@@ -0,0 +1,82 @@
+/* eg_cosine_batch_strategy.h — internal Strategy interface, NOT for call
+ * sites (they use eg_cosine_batch.h). Only eg_cosine_batch.c's factory and
+ * the concrete strategy implementation files include this.
+ *
+ * Each concrete strategy exposes exactly one getter returning a pointer to a
+ * static, immutable EgCosineBatchStrategy vtable. Which getters actually
+ * exist as linkable symbols is a BUILD-TIME concern (decided by
+ * build_vindex_bench.sh / the engram daemon's own build, via which .c/.m
+ * files get compiled per platform) gated by the EG_HAVE_STRATEGY_* macros
+ * below — the factory in eg_cosine_batch.c is the ONLY place that branches
+ * on those macros. Call sites never see them; that's the whole point of the
+ * Adapter in eg_cosine_batch.h.
+ *
+ * Three concrete strategies exist:
+ *   eg_cosine_batch_strategy_ggml()        — ggml + dynamically-loaded Metal
+ *                                             backend plugin. Darwin only in
+ *                                             this build; the default
+ *                                             preferred strategy wherever
+ *                                             available. EG_HAVE_STRATEGY_GGML.
+ *   eg_cosine_batch_strategy_metal_hand()  — the original hand-rolled Metal
+ *                                             compute shader from PR #114
+ *                                             (eg_cosine_batch.metal),
+ *                                             preserved verbatim as a
+ *                                             selectable fallback strategy,
+ *                                             not deleted. Darwin only.
+ *                                             EG_HAVE_STRATEGY_METAL_HAND.
+ *   eg_cosine_batch_strategy_cpu()         — universal always-false
+ *                                             fallback. Always compiled, on
+ *                                             every platform; this is what a
+ *                                             non-Darwin build links
+ *                                             exclusively (matching PR #114's
+ *                                             eg_metal_cosine_stub.c), and
+ *                                             what any platform falls back
+ *                                             to when no real strategy is
+ *                                             available at runtime.
+ */
+#ifndef EG_COSINE_BATCH_STRATEGY_H
+#define EG_COSINE_BATCH_STRATEGY_H
+
+#include <stdint.h>
+#include <stdbool.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct EgCosineBatchStrategy {
+    /* Stable, short, lowercase-hyphenated identifier — what
+     * eg_cosine_batch_strategy_name() surfaces. Never NULL. */
+    const char* name;
+
+    /* Cheap after the first call (lazy init, cached internally). Must never
+     * throw/crash/hang — mirrors eg_cosine_batch_available()'s contract. */
+    bool (*available)(void);
+
+    /* Same shape/contract as eg_cosine_batch() in eg_cosine_batch.h. */
+    bool (*batch)(const float* query, int32_t qdim,
+                  const float* const* node_ptrs, const int32_t* node_dims,
+                  int32_t n, double* out_scores);
+
+    /* Same shape/contract as eg_cosine_batch_multi() in eg_cosine_batch.h. */
+    bool (*batch_multi)(const float* queries, int32_t qdim, int32_t nq,
+                         const float* const* node_ptrs, const int32_t* node_dims,
+                         int32_t n, double* out_scores);
+} EgCosineBatchStrategy;
+
+#ifdef EG_HAVE_STRATEGY_GGML
+const EgCosineBatchStrategy* eg_cosine_batch_strategy_ggml(void);
+#endif
+
+#ifdef EG_HAVE_STRATEGY_METAL_HAND
+const EgCosineBatchStrategy* eg_cosine_batch_strategy_metal_hand(void);
+#endif
+
+/* Always declared/linked, on every platform/build. */
+const EgCosineBatchStrategy* eg_cosine_batch_strategy_cpu(void);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* EG_COSINE_BATCH_STRATEGY_H */

lang/runtime/eg_cosine_batch_strategy_cpu.c

@@ -0,0 +1,45 @@
+/* eg_cosine_batch_strategy_cpu.c — plain-C, zero-dependency universal
+ * fallback strategy. Always returns false / unavailable. Direct descendant
+ * of PR #114's eg_metal_cosine_stub.c, generalized from "the Metal stub" to
+ * "the strategy vtable's universal fallback entry" now that multiple real
+ * strategies can exist.
+ *
+ * Always compiled, on every platform. On Darwin builds it is the last-resort
+ * strategy the factory falls back to when neither ggml nor the hand-rolled
+ * Metal strategy is available at runtime (no device, compile failure, ...).
+ * On non-Darwin builds it is the ONLY strategy compiled in at all — no
+ * Objective-C, no Metal frameworks, no ggml/Metal backend plugin — so
+ * eg_cosine_batch()/eg_cosine_batch_multi() always return false there and
+ * every call site's existing CPU fallback runs unconditionally, exactly as
+ * before this PR.
+ */
+#include "eg_cosine_batch_strategy.h"
+
+static bool cpu_available(void) {
+    return false;
+}
+
+static bool cpu_batch(const float* query, int32_t qdim,
+                       const float* const* node_ptrs, const int32_t* node_dims,
+                       int32_t n, double* out_scores) {
+    (void)query; (void)qdim; (void)node_ptrs; (void)node_dims; (void)n; (void)out_scores;
+    return false;
+}
+
+static bool cpu_batch_multi(const float* queries, int32_t qdim, int32_t nq,
+                             const float* const* node_ptrs, const int32_t* node_dims,
+                             int32_t n, double* out_scores) {
+    (void)queries; (void)qdim; (void)nq; (void)node_ptrs; (void)node_dims; (void)n; (void)out_scores;
+    return false;
+}
+
+static const EgCosineBatchStrategy g_cpu_strategy = {
+    .name = "cpu-fallback",
+    .available = cpu_available,
+    .batch = cpu_batch,
+    .batch_multi = cpu_batch_multi,
+};
+
+const EgCosineBatchStrategy* eg_cosine_batch_strategy_cpu(void) {
+    return &g_cpu_strategy;
+}

lang/runtime/eg_cosine_batch_strategy_ggml.c

@@ -0,0 +1,515 @@
+/* eg_cosine_batch_strategy_ggml.c — the GGML Strategy, and the preferred
+ * default whenever it is available (see the factory's selection order in
+ * eg_cosine_batch.c).
+ *
+ * WHY: directive from Will Anderson — stop hand-rolling GPU kernels, use a
+ * real, proven, permissively-licensed library instead. ggml (the compute
+ * library underneath llama.cpp, MIT licensed) is already installed on this
+ * machine as a standalone Homebrew package (`brew info ggml`), independent
+ * of llama.cpp itself. This file is a genuinely bounded COMPUTE UTILITY —
+ * batch cosine-similarity math — analogous to a VBD Accessor calling out to
+ * infrastructure. It is explicitly NOT the engram's reasoning/persistence
+ * core; using ggml here does not cross the "own the core" line, because
+ * batch cosine math is infrastructure, not the graph traversal / activation
+ * spreading / "thinking" that IS the core and stays 100% own-code.
+ *
+ * ── The real API shape (verified against the installed headers + a
+ *    standalone probe program, not assumed from memory of other tensor
+ *    libraries) ──────────────────────────────────────────────────────────
+ *
+ * ggml ships its CPU and Metal implementations as DYNAMICALLY LOADED PLUGIN
+ * .so files (confirmed by nm: `ggml_backend_metal_init` is NOT an exported
+ * symbol of libggml.dylib/libggml-base.dylib — it exists ONLY inside
+ * libggml-metal.so under $(brew --prefix ggml)/libexec/). You cannot link
+ * `-lggml-metal`; you must go through ggml's backend REGISTRY:
+ *
+ *   1. ggml_backend_load_all_from_path(dir) — dlopen()s every backend plugin
+ *      .so found in `dir` and registers its device(s). We point this at
+ *      $(brew --prefix ggml)/libexec (resolved once, at build+init time; see
+ *      eg_ggml_backend_dir() below) rather than relying on
+ *      ggml_backend_load_all()'s own default search heuristics, which are
+ *      tuned for an installed llama.cpp-style app bundle layout, not an
+ *      arbitrary `cc`-built binary invoked from an arbitrary cwd — the exact
+ *      same "must not silently fall back to CPU for reasons that have
+ *      nothing to do with GPU availability" concern PR #114's hand-rolled
+ *      bridge already documented for its own embedded-shader-source choice.
+ *   2. ggml_backend_dev_by_type(GGML_BACKEND_DEVICE_TYPE_GPU) — find the
+ *      registered Metal device.
+ *   3. ggml_backend_dev_init(dev, NULL) — get a live ggml_backend_t.
+ *   4. Build a tiny ggml_context (no_alloc=true; it holds only tensor
+ *      metadata, not data), declare 2D F32 tensors, ggml_mul_mat(nodes,
+ *      query) — ggml's documented convention: A is [k cols, n rows], B is
+ *      [k cols, m rows] (transposed internally), result is [n cols, m rows]
+ *      — i.e. mul_mat(node_matrix[dim,n], query_matrix[dim,nq]) yields
+ *      out[n,nq] where out[j*n+i] = dot(node_i, query_j). A row-major
+ *      (dim,n) node matrix and a row-major (dim,nq) query matrix is EXACTLY
+ *      the packed layout the hand-rolled Metal kernel already used — one
+ *      matmul replaces the whole per-row dot-product loop.
+ *   5. ggml_backend_alloc_ctx_tensors(ctx, backend) to actually allocate
+ *      device buffers for those tensors, ggml_backend_tensor_set() to upload,
+ *      ggml_backend_graph_compute() to run, ggml_backend_tensor_get() to
+ *      read back.
+ *
+ * This exact sequence was verified end-to-end in a standalone probe (build
+ * it yourself: see the PR description) against a plain-C CPU dot product —
+ * bit-for-bit correct within float rounding. Real numbers against the
+ * el_runtime.c CPU oracle are reported in the PR body via vindex_bench.
+ *
+ * ggml_mul_mat only computes the raw dot products — it has no notion of
+ * "cosine" or of this codebase's -2.0 dim-mismatch/null/zero-norm sentinel.
+ * Per the adapter's directive: gather only VALID, uniform-dim rows into the
+ * packed matrix sent to the GPU (skipping null/mismatched rows entirely,
+ * rather than the hand-rolled kernel's zero-pad-and-sentinel-in-shader
+ * approach), then scatter -2.0 back for every row that was excluded — same
+ * gather/scatter contract eg_cosine_batch.h documents. Norms (||node||,
+ * ||query||) are computed on the CPU host in the same pass that already
+ * touches every element to gather/convert — essentially free — using the
+ * same 4-way-partial-sum accumulation the hand-rolled kernel and the CPU
+ * oracle both use, so the float32 error profile stays comparable across all
+ * three strategies. Only the O(n*dim*nq) dot-product matmul — the actual
+ * expensive part — is offloaded to the GPU.
+ *
+ * ── Precision: the ne11<=8 chunking, and why it is not optional ──────────
+ *
+ * The claim in the first version of this file — "ggml_mul_mat on F32 x F32
+ * inputs computes in F32 on the Metal backend" — is WRONG, and the 0.9933
+ * id-recall it shipped with (vs the hand-rolled kernel's 0.9997) was the
+ * symptom. ggml-metal has two F32xF32 matmul kernels and picks between them
+ * purely on ne11 (the number of B rows == our query count):
+ *
+ *   ne11 <= 8  -> kernel_mul_mv_ext_f32_f32_* / kernel_mul_mv_f32_f32_*
+ *                 templated <float, float> — genuine F32 accumulation.
+ *   ne11 >  8  -> kernel_mul_mm_f32_f32, which is templated
+ *                 <half, half4x4, simdgroup_half8x8, half, half2x4,
+ *                  simdgroup_half8x8, ...> — i.e. BOTH operands are narrowed
+ *                 to F16 and accumulated in simdgroup_half8x8 tiles, even
+ *                 though the tensors are GGML_TYPE_F32 on both sides.
+ *
+ * (Read it yourself, no guessing — the kernel templates are literal strings
+ * in the shipped plugin:
+ *    strings $(brew --prefix ggml)/libexec/libggml-metal.so \
+ *      | grep -E 'host_name\("kernel_mul_m[mv]_f32_f32'
+ *  and the runtime pick is visible with GGML_METAL_DEBUG-style logging as
+ *  "compiling pipeline: base = 'kernel_mul_mm_f32_f32'".)
+ *
+ * The previous code issued ONE ggml_mul_mat with ne11 = nq (300 in the
+ * benchmark), landing squarely on the F16 mul_mm path. Measured on this
+ * machine (M4 Pro), n=13415 x dim=768 x nq=300, against a CPU double-
+ * accumulated oracle:
+ *
+ *   ne11=300 (one mul_mat, the old code) : mean |Δdot| = 1.038e-05
+ *   ne11=8   (chunked, this code)        : mean |Δdot| = 3.863e-09
+ *
+ * — a ~2700x reduction in dot-product error, which is exactly the gap that
+ * showed up as 0.9933-vs-0.9997 recall.
+ *
+ * ggml_mul_mat_set_prec(t, GGML_PREC_F32) does NOT fix this. It was tried:
+ * the error was bit-identical with and without it (1.038e-05 either way),
+ * because ggml-metal only consults the prec flag on paths that have an F32
+ * variant to switch to, and there is no F32-accumulating mul_mm kernel in
+ * this build to select. The ONLY lever from outside ggml is ne11.
+ *
+ * So: instead of one mul_mat with ne11=nq, we emit ceil(nq/8) mul_mats, each
+ * over an ne11<=8 ggml_view_2d slice of the same query tensor, all into ONE
+ * graph and ONE ggml_backend_graph_compute. The node matrix is still uploaded
+ * exactly once and still read by the GPU as one shared operand — the whole
+ * point of batch_multi is preserved.
+ *
+ * The cost is real, and stated rather than buried. Timing the whole
+ * batch_multi() call (gather + norms + upload + GPU + scatter) on the real
+ * shape, median of 15 reps after a discarded warm-up, three separate runs:
+ *
+ *   unchunked (old, F16 mm) : 13.19 / 13.35 / 14.42 ms   -> ~0.044 ms/query
+ *   chunked   (this code)   : 19.92 / 20.08 / 20.23 ms   -> ~0.067 ms/query
+ *   hand-rolled Metal       : 17.74 / 17.88 / 17.99 ms   -> ~0.060 ms/query
+ *
+ * So correctness here costs about +6.7ms per 300-query batch (~1.5x on this
+ * call), and leaves us ~12% behind the hand-rolled kernel instead of ~35%
+ * ahead of it. That is not free and should not be sold as free. The reason it
+ * cannot be recovered inside ggml: an fp32 matmul on Metal has to re-stream
+ * the whole node matrix once per <=8 queries (38 dispatches x ~41MB here),
+ * where the F16 mul_mm kernel tiles it in threadgroup memory and reads it far
+ * fewer times. ggml's Metal backend ships no fp32 TILED matmul, so on this
+ * backend "fast" and "fp32" are genuinely exclusive — the hand-rolled kernel
+ * escapes the choice only because it is an fp32 kernel written for this one
+ * shape. Trading precision back for speed is a one-line env change; trading
+ * the other way was not available before this commit at all.
+ *
+ * 8 is not a magic number we invented — it is ggml-metal's own mul_mm
+ * threshold, measured by sweeping ne11 and watching both the error and which
+ * pipeline ggml compiles (9 flips to mul_mm and the error jumps back to
+ * 1.0e-05 in the same step). EL_GGML_MULMAT_CHUNK overrides it: raise it to
+ * trade this precision back for throughput, or set it >= nq to reproduce the
+ * old single-mul_mat behaviour exactly. If a future ggml moves the threshold,
+ * the worst case is that we silently land back on mul_mm — the same accuracy
+ * we shipped before, never a correctness break.
+ *
+ * ── Cold start: what is and is not ours to fix ───────────────────────────
+ *
+ * The ~7.8s first-call cost reported for the first version of this file is
+ * NOT this file re-initialising per call (init is, and always was, cached
+ * behind g_init_attempted below). It is Apple's Metal shader cache missing
+ * on ggml's embedded metallib — ggml-metal ships ~650 kernels in one
+ * __ggml_metallib section, and the first newLibraryWithData of it on a given
+ * machine costs seconds ("ggml_metal_library_init: loaded in 7.670 sec")
+ * while the driver populates ~/…/C/com.apple.metal/. That cache is keyed on
+ * the library, not on our binary, and is shared across processes: the very
+ * next run of a DIFFERENT binary linking the same ggml reports
+ * "loaded in 0.009 sec". So it is a once-per-machine, per-ggml-version cost,
+ * not a per-process one, and nothing this file does can avoid it — the
+ * hand-rolled strategy escapes it only because its shader is two small
+ * kernels instead of six hundred.
+ *
+ * The residual warm init IS ours to look at, and the answer there is "there
+ * was nothing much to win": ggml_backend_load_all_from_path() dlopens every
+ * plugin in the directory (three CPU micro-arch variants + BLAS + Metal) when
+ * we only ever use Metal, so we now load the single Metal plugin instead —
+ * but measured warm that is 44.7-52.4ms against 46.9-58.9ms, i.e. the same
+ * number inside noise, because libggml-metal.so's own init dominates. Warm
+ * ggml init lands at 44-53ms, against 36-117ms for the hand-rolled strategy's
+ * device+pipeline setup. Cold start was never the real defect here; precision
+ * was.
+ */
+#include "eg_cosine_batch_strategy.h"
+
+#include <ggml.h>
+#include <ggml-backend.h>
+#include <ggml-alloc.h>
+
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+/* ── lazy, one-time backend init, cached ─────────────────────────────────── */
+static bool g_init_attempted = false;
+static bool g_init_ok = false;
+static ggml_backend_t g_backend = NULL;
+
+/* Where to look for the dynamically-loaded backend plugin .so files.
+ * EL_GGML_BACKEND_PATH overrides for non-standard installs; otherwise we try
+ * the Homebrew opt-prefix symlink (stable across ggml point-version bumps —
+ * $(brew --prefix ggml)/libexec — confirmed to exist and contain
+ * libggml-metal.so / libggml-cpu-*.so / libggml-blas.so on this machine),
+ * falling back to ggml's own default search (ggml_backend_load_all()) in
+ * case a different install layout (e.g. a from-source build with a
+ * standard-prefix install) makes that succeed instead. */
+static const char* eg_ggml_backend_dir(void) {
+    const char* s = getenv("EL_GGML_BACKEND_PATH");
+    if (s && *s) return s;
+    return "/opt/homebrew/opt/ggml/libexec";
+}
+
+/* "<dir>/libggml-metal.so" in a static buffer. Only ever called once, from
+ * eg_ggml_ensure_init(), before any thread could race it. */
+static const char* eg_ggml_metal_plugin_path(const char* dir) {
+    static char buf[1024];
+    snprintf(buf, sizeof buf, "%s/libggml-metal.so", dir);
+    return buf;
+}
+
+/* Largest ne11 (query-batch rows per ggml_mul_mat) that keeps ggml-metal on
+ * its F32 mul_mv kernels instead of the F16-accumulating mul_mm kernel — see
+ * the precision discussion in this file's header. EL_GGML_MULMAT_CHUNK
+ * overrides; a value <= 0 means "use the default". */
+#define EG_GGML_MULMAT_CHUNK_DEFAULT 8
+
+static int32_t eg_ggml_mulmat_chunk(void) {
+    static bool resolved = false;
+    static int32_t chunk = EG_GGML_MULMAT_CHUNK_DEFAULT;
+    if (!resolved) {
+        resolved = true;
+        const char* s = getenv("EL_GGML_MULMAT_CHUNK");
+        if (s && *s) {
+            long v = strtol(s, NULL, 10);
+            if (v > 0 && v <= INT32_MAX) chunk = (int32_t)v;
+        }
+    }
+    return chunk;
+}
+
+/* Which ggml device this strategy computes on. GPU (Metal) is the default
+ * because offloading is the architectural point — the engram's own graph
+ * traversal and activation spreading are CPU work, and a "GPU" strategy that
+ * quietly saturates the CPU steals from them.
+ *
+ * ACCEL (ggml's BLAS/Accelerate plugin) is reachable here mainly as a
+ * portability fallback and a diagnostic, and it is documented as MEASURED AND
+ * REJECTED rather than as a recommendation. In an isolated probe that timed
+ * only ggml_backend_graph_compute, BLAS looked excellent — 3.4-4.0ms for the
+ * 300-query batch at mean |Δdot| 1.5e-08, i.e. as fast as the old F16 path and
+ * far more accurate. End to end on the real store through vindex_bench it does
+ * not hold up: 0.191 ms/query at id-recall 0.9973, against 0.125-0.142 ms/query
+ * at 0.9987 for the Metal default. It is dominated on BOTH axes, because the
+ * isolated probe was not competing with the rest of the batch for the same CPU
+ * cores and the real call path is. Kept because a machine with no usable Metal
+ * device still wants a working ggml strategy — not because it is faster. */
+static enum ggml_backend_dev_type eg_ggml_device_type(void) {
+    const char* s = getenv("EL_GGML_DEVICE");
+    if (s && *s) {
+        if (strcmp(s, "accel") == 0) return GGML_BACKEND_DEVICE_TYPE_ACCEL;
+        if (strcmp(s, "cpu")   == 0) return GGML_BACKEND_DEVICE_TYPE_CPU;
+    }
+    return GGML_BACKEND_DEVICE_TYPE_GPU;
+}
+
+static bool eg_ggml_ensure_init(void) {
+    if (g_init_attempted) return g_init_ok;
+    g_init_attempted = true;
+
+    const char* dir = eg_ggml_backend_dir();
+    const enum ggml_backend_dev_type want = eg_ggml_device_type();
+
+    /* Metal is the only backend this strategy uses by default, so load just
+     * that one plugin rather than dlopening the whole directory (three CPU
+     * micro-arch variants + BLAS + Metal here).
+     *
+     * Be honest about what this buys: almost nothing in wall time. Measured
+     * warm, three runs each — load-everything 58.9/46.9/55.1ms, Metal-only
+     * 52.4/51.2/44.7ms. The cost is dominated by dlopening and initialising
+     * libggml-metal.so itself, not by the four plugins we skip, so the two
+     * overlap inside noise. It is kept because registering four device types
+     * we will never dispatch to is untidy and makes ggml_backend_dev_by_type
+     * ambiguous, not because it is a speedup — do not cite it as one.
+     *
+     * ggml_backend_load() returns NULL for a missing or unloadable path,
+     * which simply falls through to the broader searches below; it is never
+     * fatal. Any non-default device needs the full directory scan to find
+     * its plugin, so skip the fast path there. */
+    if (want == GGML_BACKEND_DEVICE_TYPE_GPU)
+        ggml_backend_load(eg_ggml_metal_plugin_path(dir));
+
+    ggml_backend_dev_t dev = ggml_backend_dev_by_type(want);
+    if (!dev) {
+        /* Non-standard layout, a ggml built with a differently-named Metal
+         * plugin, or a non-default device: dlopen every plugin in `dir`. */
+        ggml_backend_load_all_from_path(dir);
+        dev = ggml_backend_dev_by_type(want);
+    }
+    if (!dev) {
+        /* Fall back to ggml's own default search heuristics only if the
+         * explicit path above found nothing — avoids double-registering the
+         * same plugins (ggml does not dedupe two different paths that
+         * happen to resolve to the same files, e.g. our stable opt-prefix
+         * symlink vs. its own Cellar-relative guess) in the common case
+         * where the explicit path already worked. */
+        ggml_backend_load_all();
+        dev = ggml_backend_dev_by_type(want);
+    }
+    if (!dev) return false;
+
+    ggml_backend_t backend = ggml_backend_dev_init(dev, NULL);
+    if (!backend) return false;
+
+    g_backend = backend;
+    g_init_ok = true;
+    return true;
+}
+
+static bool ggml_strategy_available(void) {
+    return eg_ggml_ensure_init();
+}
+
+/* ── shared core: gather valid rows + norms, matmul, scatter ────────────── */
+
+/* 4-way partial-sum squared-norm accumulation over `dim` floats — same shape
+ * as eg_cosine_batch.metal's per-thread accumulation and vindex_bench.c's
+ * CPU brute_topk unroll, kept consistent on purpose so the float32 error
+ * profile is comparable across all three strategies. */
+static float eg_norm_sq_f32(const float* v, int32_t dim) {
+    float s0 = 0, s1 = 0, s2 = 0, s3 = 0;
+    int32_t d = 0, dim4 = dim & ~3;
+    for (; d < dim4; d += 4) {
+        s0 += v[d] * v[d]; s1 += v[d+1] * v[d+1];
+        s2 += v[d+2] * v[d+2]; s3 += v[d+3] * v[d+3];
+    }
+    float s = (s0 + s1) + (s2 + s3);
+    for (; d < dim; d++) s += v[d] * v[d];
+    return s;
+}
+
+/* Runs one ggml_mul_mat(node_matrix[dim,n_valid], query_matrix[dim,nq]) and
+ * combines it with CPU-computed norms into cosine scores, scattering into
+ * out_scores at ORIGINAL (ungathered) indices. out_scores must already be
+ * fully sized for n*nq (or n for the single-query case, nq=1) — every entry
+ * gets written (valid rows get a real cosine, invalid rows get -2.0), so
+ * this never leaves a partial result. Returns false only on a genuine
+ * failure (alloc, compute) at which point out_scores is left as whatever a
+ * caller-supplied scratch buffer already contained — callers here always
+ * pass a fresh buffer they discard on false, matching the adapter contract
+ * of "on failure, out_scores is treated as untouched" from the caller's
+ * point of view. */
+static bool eg_ggml_run(const float* queries, int32_t qdim, int32_t nq,
+                         const float* const* node_ptrs, const int32_t* node_dims,
+                         int32_t n, double* out_scores) {
+    if (!queries || qdim <= 0 || nq <= 0 || !node_ptrs || !node_dims || n <= 0 || !out_scores)
+        return false;
+    if (!eg_ggml_ensure_init()) return false;
+
+    /* Pass 1 (CPU): gather valid rows (non-NULL ptr, dim == qdim) into a
+     * packed (dim, n_valid) row-major matrix, remembering the original index
+     * of each packed row, and compute each valid row's squared norm in the
+     * same pass. Rows excluded here get -2.0 scattered for every query
+     * below without ever touching the GPU. */
+    int32_t* valid_orig = (int32_t*)malloc((size_t)n * sizeof(int32_t));
+    float* node_norm_sq  = (float*)malloc((size_t)n * sizeof(float)); /* indexed by packed position */
+    float* node_matrix   = NULL;
+    if (!valid_orig || !node_norm_sq) { free(valid_orig); free(node_norm_sq); return false; }
+
+    int32_t n_valid = 0;
+    for (int32_t i = 0; i < n; i++) {
+        if (node_ptrs[i] && node_dims[i] == qdim) n_valid++;
+    }
+
+    if (n_valid > 0) {
+        node_matrix = (float*)malloc((size_t)n_valid * (size_t)qdim * sizeof(float));
+        if (!node_matrix) { free(valid_orig); free(node_norm_sq); return false; }
+        int32_t w = 0;
+        for (int32_t i = 0; i < n; i++) {
+            if (!node_ptrs[i] || node_dims[i] != qdim) continue;
+            memcpy(node_matrix + (size_t)w * qdim, node_ptrs[i], (size_t)qdim * sizeof(float));
+            node_norm_sq[w] = eg_norm_sq_f32(node_ptrs[i], qdim);
+            valid_orig[w] = i;
+            w++;
+        }
+    }
+
+    /* Query norms — nq is typically small (1 or the size of one batch of
+     * comparison queries), so this loop is cheap regardless. */
+    float* q_norm_sq = (float*)malloc((size_t)nq * sizeof(float));
+    if (!q_norm_sq) { free(valid_orig); free(node_norm_sq); free(node_matrix); return false; }
+    for (int32_t j = 0; j < nq; j++) q_norm_sq[j] = eg_norm_sq_f32(queries + (size_t)j * qdim, qdim);
+
+    /* Nothing valid to compare against: every output is -2.0. Still a fully
+     * and correctly populated result — no GPU dispatch was needed to know
+     * that. */
+    if (n_valid == 0) {
+        for (size_t k = 0; k < (size_t)n * (size_t)nq; k++) out_scores[k] = -2.0;
+        free(valid_orig); free(node_norm_sq); free(node_matrix); free(q_norm_sq);
+        return true;
+    }
+
+    /* Pass 2 (GPU via ggml): dot[j*n_valid + i] = dot(node_i, query_j),
+     * computed as ceil(nq/chunk) separate ggml_mul_mat ops over ne11<=chunk
+     * ggml_view_2d slices of ONE query tensor, all expanded into ONE graph
+     * and run by ONE ggml_backend_graph_compute. Chunking is what keeps
+     * ggml-metal on its F32 mul_mv kernels rather than the F16-accumulating
+     * mul_mm kernel (see this file's header); sharing one graph and one
+     * t_nodes tensor is what keeps the node matrix uploaded exactly once,
+     * which is the entire reason batch_multi exists. */
+    const int32_t chunk = eg_ggml_mulmat_chunk();
+    const int32_t ngroups = (nq + chunk - 1) / chunk;
+
+    /* Tensors held by the context: t_nodes, t_query, plus one view and one
+     * mul_mat result per group. The graph holds at most one node per view and
+     * one per mul_mat. Slack on both so a ggml that bookkeeps slightly
+     * differently cannot silently overflow the arena. */
+    const size_t n_tensors  = (size_t)2 * (size_t)ngroups + 8;
+    const size_t graph_size = (size_t)2 * (size_t)ngroups + 16;
+    struct ggml_init_params gp = {
+        .mem_size   = ggml_tensor_overhead() * n_tensors
+                    + ggml_graph_overhead_custom(graph_size, false),
+        .mem_buffer = NULL,
+        .no_alloc   = true,
+    };
+    struct ggml_context* ctx = ggml_init(gp);
+    if (!ctx) { free(valid_orig); free(node_norm_sq); free(node_matrix); free(q_norm_sq); return false; }
+
+    struct ggml_tensor** t_dots = (struct ggml_tensor**)malloc((size_t)ngroups * sizeof(*t_dots));
+    if (!t_dots) { ggml_free(ctx); free(valid_orig); free(node_norm_sq); free(node_matrix); free(q_norm_sq); return false; }
+
+    struct ggml_tensor* t_nodes = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, qdim, n_valid);
+    struct ggml_tensor* t_query = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, qdim, nq);
+    struct ggml_cgraph*  gf     = t_nodes && t_query
+                                ? ggml_new_graph_custom(ctx, graph_size, false) : NULL;
+    if (!gf) { free(t_dots); ggml_free(ctx); free(valid_orig); free(node_norm_sq); free(node_matrix); free(q_norm_sq); return false; }
+
+    bool built = true;
+    for (int32_t g = 0; g < ngroups; g++) {
+        const int32_t start = g * chunk;
+        const int32_t count = (start + chunk <= nq) ? chunk : (nq - start);
+        struct ggml_tensor* t_qv = ggml_view_2d(ctx, t_query, qdim, count,
+                                                t_query->nb[1],
+                                                (size_t)start * t_query->nb[1]);
+        t_dots[g] = t_qv ? ggml_mul_mat(ctx, t_nodes, t_qv) : NULL;
+        if (!t_dots[g]) { built = false; break; }
+        ggml_build_forward_expand(gf, t_dots[g]);
+    }
+    if (!built) { free(t_dots); ggml_free(ctx); free(valid_orig); free(node_norm_sq); free(node_matrix); free(q_norm_sq); return false; }
+
+    struct ggml_backend_buffer* buf = ggml_backend_alloc_ctx_tensors(ctx, g_backend);
+    if (!buf) { free(t_dots); ggml_free(ctx); free(valid_orig); free(node_norm_sq); free(node_matrix); free(q_norm_sq); return false; }
+
+    ggml_backend_tensor_set(t_nodes, node_matrix, 0, (size_t)n_valid * qdim * sizeof(float));
+    ggml_backend_tensor_set(t_query, queries, 0, (size_t)nq * qdim * sizeof(float));
+    free(node_matrix); /* uploaded; the packed CPU copy is no longer needed */
+
+    enum ggml_status st = ggml_backend_graph_compute(g_backend, gf);
+    if (st != GGML_STATUS_SUCCESS) {
+        free(t_dots); ggml_backend_buffer_free(buf); ggml_free(ctx);
+        free(valid_orig); free(node_norm_sq); free(q_norm_sq);
+        return false;
+    }
+
+    /* Each group's result is [n_valid, count] contiguous, so reading group g
+     * into dot + start*n_valid reconstructs exactly the same flat
+     * dot[j*n_valid + w] layout a single ne11=nq mul_mat would have produced —
+     * Pass 3 below is unchanged by the chunking. */
+    float* dot = (float*)malloc((size_t)n_valid * (size_t)nq * sizeof(float));
+    if (!dot) { free(t_dots); ggml_backend_buffer_free(buf); ggml_free(ctx); free(valid_orig); free(node_norm_sq); free(q_norm_sq); return false; }
+    for (int32_t g = 0; g < ngroups; g++) {
+        const int32_t start = g * chunk;
+        const int32_t count = (start + chunk <= nq) ? chunk : (nq - start);
+        ggml_backend_tensor_get(t_dots[g], dot + (size_t)start * n_valid, 0,
+                                (size_t)count * (size_t)n_valid * sizeof(float));
+    }
+    free(t_dots);
+
+    /* Pass 3 (CPU): combine dot/(||a||*||b||) per (query,node) pair, scatter
+     * into out_scores at ORIGINAL node indices; every excluded row gets
+     * -2.0 for every query. out_scores is fully populated either way. */
+    for (int32_t j = 0; j < nq; j++) {
+        double* orow = out_scores + (size_t)j * n;
+        for (int32_t i = 0; i < n; i++) orow[i] = -2.0; /* default: excluded */
+        for (int32_t w = 0; w < n_valid; w++) {
+            float na = node_norm_sq[w], nb = q_norm_sq[j];
+            int32_t oi = valid_orig[w];
+            if (na <= 0.0f || nb <= 0.0f) { orow[oi] = -2.0; continue; }
+            float d = dot[(size_t)j * n_valid + w];
+            orow[oi] = (double)(d / sqrtf(na * nb));
+        }
+    }
+
+    free(dot);
+    ggml_backend_buffer_free(buf);
+    ggml_free(ctx);
+    free(valid_orig); free(node_norm_sq); free(q_norm_sq);
+    return true;
+}
+
+static bool ggml_strategy_batch(const float* query, int32_t qdim,
+                                 const float* const* node_ptrs, const int32_t* node_dims,
+                                 int32_t n, double* out_scores) {
+    if (!query || qdim <= 0 || !node_ptrs || !node_dims || n <= 0 || !out_scores) return false;
+    /* out_scores here is n doubles (nq=1); eg_ggml_run writes n*nq = n of
+     * them, laid out identically to the single-query contract. */
+    return eg_ggml_run(query, qdim, 1, node_ptrs, node_dims, n, out_scores);
+}
+
+static bool ggml_strategy_batch_multi(const float* queries, int32_t qdim, int32_t nq,
+                                       const float* const* node_ptrs, const int32_t* node_dims,
+                                       int32_t n, double* out_scores) {
+    return eg_ggml_run(queries, qdim, nq, node_ptrs, node_dims, n, out_scores);
+}
+
+static const EgCosineBatchStrategy g_ggml_strategy = {
+    .name = "ggml",
+    .available = ggml_strategy_available,
+    .batch = ggml_strategy_batch,
+    .batch_multi = ggml_strategy_batch_multi,
+};
+
+const EgCosineBatchStrategy* eg_cosine_batch_strategy_ggml(void) {
+    return &g_ggml_strategy;
+}

lang/runtime/eg_cosine_batch_strategy_metal_hand.m

@@ -0,0 +1,358 @@
+/* eg_cosine_batch_strategy_metal_hand.m — the HAND-ROLLED-METAL Strategy.
+ *
+ * This is PR #114's original Objective-C bridge (formerly eg_metal_cosine.m)
+ * exposing the hand-written Metal compute shader (eg_cosine_batch.metal) as
+ * one concrete EgCosineBatchStrategy. It is preserved here almost verbatim —
+ * real, carefully verified work, not discarded — now living behind the
+ * Adapter/Strategy/Factory restructuring (see eg_cosine_batch.h and
+ * eg_cosine_batch_strategy.h) alongside the new ggml-Metal strategy
+ * (eg_cosine_batch_strategy_ggml.c) and the universal CPU fallback
+ * (eg_cosine_batch_strategy_cpu.c). The factory in eg_cosine_batch.c prefers
+ * ggml by default when both are available; this strategy remains selectable
+ * via EL_COSINE_BATCH_STRATEGY=metal, and is what the factory falls back to
+ * if ggml's backend plugin fails to load/init for any reason.
+ *
+ * Apple-only (Metal has no other platform). This file is excluded from the
+ * build entirely on non-Darwin — see build_vindex_bench.sh, which only
+ * compiles/links this file and defines EG_HAVE_STRATEGY_METAL_HAND when
+ * `uname` is Darwin. On Linux the factory never sees this strategy at all —
+ * callers must always be prepared for the "no real strategy available"
+ * fallback via the CPU strategy, which is also exactly what happens here on
+ * Apple hardware with no usable GPU.
+ *
+ * Design (unchanged from PR #114):
+ *   - Device/queue/pipeline are created lazily, once, and cached in static
+ *     globals — every call after the first only allocates buffers + submits.
+ *   - The Metal shader source is embedded as a C string literal (kMetalSrc
+ *     below) rather than loaded from a file at runtime or shipped as a
+ *     precompiled .metallib. Chosen over newLibraryWithFile: /a .metallib
+ *     because the engram binary can be invoked from an arbitrary working
+ *     directory (launchd job, nsbx sandbox, CI) and a file-path shader would
+ *     be one relocation away from silently falling back to CPU for reasons
+ *     that have nothing to do with Metal availability. Embedding costs one
+ *     runtime shader compile (~tens of ms) on first use, amortized over the
+ *     process lifetime, in exchange for a genuinely self-contained binary.
+ *     Source of truth for review/tooling is eg_cosine_batch.metal — this
+ *     string MUST be kept byte-identical to that file (a comment marks both
+ *     ends of the copy).
+ *   - Buffers use MTLResourceStorageModeShared: on Apple Silicon's unified
+ *     memory, CPU and GPU read the same physical pages, so filling a buffer
+ *     is a plain memcpy and there is no separate "upload" step.
+ *   - ANY failure at ANY step (no device, pipeline compile error, buffer
+ *     allocation failure, bad args) returns false and leaves out_scores
+ *     untouched. This function is called from the request-handling hot path
+ *     of a long-lived daemon — it must never throw, crash, or hang it.
+ */
+#import <Foundation/Foundation.h>
+#import <Metal/Metal.h>
+#include "eg_cosine_batch_strategy.h"
+#include <string.h>
+#include <stdlib.h>
+
+/* ── BEGIN embedded shader source (keep in sync with eg_cosine_batch.metal) ── */
+static const char* kEgCosineBatchMetalSrc =
+"#include <metal_stdlib>\n"
+"using namespace metal;\n"
+"struct EgCosineParams { uint n; uint dim; };\n"
+"kernel void eg_cosine_batch_kernel(\n"
+"    device const float*  query        [[buffer(0)]],\n"
+"    device const float*  node_matrix  [[buffer(1)]],\n"
+"    device const int*    node_dims    [[buffer(2)]],\n"
+"    constant EgCosineParams& p        [[buffer(3)]],\n"
+"    device float*        out_scores   [[buffer(4)]],\n"
+"    uint gid [[thread_position_in_grid]])\n"
+"{\n"
+"    if (gid >= p.n) return;\n"
+"    if (node_dims[gid] != int(p.dim)) { out_scores[gid] = -2.0f; return; }\n"
+"    device const float* row = node_matrix + (uint64_t)gid * (uint64_t)p.dim;\n"
+"    float dot0 = 0.0f, dot1 = 0.0f, dot2 = 0.0f, dot3 = 0.0f;\n"
+"    float na0  = 0.0f, na1  = 0.0f, na2  = 0.0f, na3  = 0.0f;\n"
+"    float nb0  = 0.0f, nb1  = 0.0f, nb2  = 0.0f, nb3  = 0.0f;\n"
+"    uint d = 0;\n"
+"    uint dim4 = p.dim & ~3u;\n"
+"    for (; d < dim4; d += 4) {\n"
+"        float a0 = row[d],   b0 = query[d];\n"
+"        float a1 = row[d+1], b1 = query[d+1];\n"
+"        float a2 = row[d+2], b2 = query[d+2];\n"
+"        float a3 = row[d+3], b3 = query[d+3];\n"
+"        dot0 += a0*b0; dot1 += a1*b1; dot2 += a2*b2; dot3 += a3*b3;\n"
+"        na0  += a0*a0; na1  += a1*a1; na2  += a2*a2; na3  += a3*a3;\n"
+"        nb0  += b0*b0; nb1  += b1*b1; nb2  += b2*b2; nb3  += b3*b3;\n"
+"    }\n"
+"    float dot = (dot0 + dot1) + (dot2 + dot3);\n"
+"    float na  = (na0  + na1)  + (na2  + na3);\n"
+"    float nb  = (nb0  + nb1)  + (nb2  + nb3);\n"
+"    for (; d < p.dim; d++) {\n"
+"        float a = row[d], b = query[d];\n"
+"        dot += a*b; na += a*a; nb += b*b;\n"
+"    }\n"
+"    if (na <= 0.0f || nb <= 0.0f) { out_scores[gid] = -2.0f; return; }\n"
+"    out_scores[gid] = dot / sqrt(na * nb);\n"
+"}\n"
+"struct EgCosineMultiParams { uint n; uint dim; uint nq; };\n"
+"kernel void eg_cosine_batch_multi_kernel(\n"
+"    device const float*  queries      [[buffer(0)]],\n"
+"    device const float*  node_matrix  [[buffer(1)]],\n"
+"    device const int*    node_dims    [[buffer(2)]],\n"
+"    constant EgCosineMultiParams& p   [[buffer(3)]],\n"
+"    device float*        out_scores   [[buffer(4)]],\n"
+"    uint2 gid [[thread_position_in_grid]])\n"
+"{\n"
+"    uint nid = gid.x, qid = gid.y;\n"
+"    if (nid >= p.n || qid >= p.nq) return;\n"
+"    uint64_t out_idx = (uint64_t)qid * (uint64_t)p.n + (uint64_t)nid;\n"
+"    if (node_dims[nid] != int(p.dim)) { out_scores[out_idx] = -2.0f; return; }\n"
+"    device const float* row   = node_matrix + (uint64_t)nid * (uint64_t)p.dim;\n"
+"    device const float* query = queries     + (uint64_t)qid * (uint64_t)p.dim;\n"
+"    float dot0 = 0.0f, dot1 = 0.0f, dot2 = 0.0f, dot3 = 0.0f;\n"
+"    float na0  = 0.0f, na1  = 0.0f, na2  = 0.0f, na3  = 0.0f;\n"
+"    float nb0  = 0.0f, nb1  = 0.0f, nb2  = 0.0f, nb3  = 0.0f;\n"
+"    uint d = 0;\n"
+"    uint dim4 = p.dim & ~3u;\n"
+"    for (; d < dim4; d += 4) {\n"
+"        float a0 = row[d],   b0 = query[d];\n"
+"        float a1 = row[d+1], b1 = query[d+1];\n"
+"        float a2 = row[d+2], b2 = query[d+2];\n"
+"        float a3 = row[d+3], b3 = query[d+3];\n"
+"        dot0 += a0*b0; dot1 += a1*b1; dot2 += a2*b2; dot3 += a3*b3;\n"
+"        na0  += a0*a0; na1  += a1*a1; na2  += a2*a2; na3  += a3*a3;\n"
+"        nb0  += b0*b0; nb1  += b1*b1; nb2  += b2*b2; nb3  += b3*b3;\n"
+"    }\n"
+"    float dot = (dot0 + dot1) + (dot2 + dot3);\n"
+"    float na  = (na0  + na1)  + (na2  + na3);\n"
+"    float nb  = (nb0  + nb1)  + (nb2  + nb3);\n"
+"    for (; d < p.dim; d++) {\n"
+"        float a = row[d], b = query[d];\n"
+"        dot += a*b; na += a*a; nb += b*b;\n"
+"    }\n"
+"    if (na <= 0.0f || nb <= 0.0f) { out_scores[out_idx] = -2.0f; return; }\n"
+"    out_scores[out_idx] = dot / sqrt(na * nb);\n"
+"}\n";
+/* ── END embedded shader source ── */
+
+typedef struct EgCosineParamsC { uint32_t n; uint32_t dim; } EgCosineParamsC;
+typedef struct EgCosineMultiParamsC { uint32_t n; uint32_t dim; uint32_t nq; } EgCosineMultiParamsC;
+
+static id<MTLDevice>              g_device   = nil;
+static id<MTLCommandQueue>        g_queue    = nil;
+static id<MTLComputePipelineState> g_pipeline = nil;       /* single-query kernel */
+static id<MTLComputePipelineState> g_pipeline_multi = nil; /* multi-query kernel */
+static bool g_init_attempted = false;
+static bool g_init_ok        = false;
+
+/* Lazy, one-time setup. Never throws — every Metal call here is the
+ * "returns nil/NSError on failure" flavor, not an exception-throwing one. */
+static bool eg_metal_ensure_init(void) {
+    if (g_init_attempted) return g_init_ok;
+    g_init_attempted = true;
+
+    @autoreleasepool {
+        id<MTLDevice> dev = MTLCreateSystemDefaultDevice();
+        if (!dev) return false;
+
+        id<MTLCommandQueue> q = [dev newCommandQueue];
+        if (!q) return false;
+
+        NSError* err = nil;
+        NSString* src = [NSString stringWithUTF8String:kEgCosineBatchMetalSrc];
+        MTLCompileOptions* opts = [MTLCompileOptions new];
+        id<MTLLibrary> lib = [dev newLibraryWithSource:src options:opts error:&err];
+        if (!lib) return false;
+
+        id<MTLFunction> fn = [lib newFunctionWithName:@"eg_cosine_batch_kernel"];
+        if (!fn) return false;
+        id<MTLComputePipelineState> pipe = [dev newComputePipelineStateWithFunction:fn error:&err];
+        if (!pipe) return false;
+
+        id<MTLFunction> fnMulti = [lib newFunctionWithName:@"eg_cosine_batch_multi_kernel"];
+        if (!fnMulti) return false;
+        id<MTLComputePipelineState> pipeMulti = [dev newComputePipelineStateWithFunction:fnMulti error:&err];
+        if (!pipeMulti) return false;
+
+        g_device = dev;
+        g_queue = q;
+        g_pipeline = pipe;
+        g_pipeline_multi = pipeMulti;
+        g_init_ok = true;
+        return true;
+    }
+}
+
+static bool mh_available(void) {
+    return eg_metal_ensure_init();
+}
+
+static bool mh_batch(const float* query, int32_t qdim,
+                      const float* const* node_ptrs,
+                      const int32_t* node_dims,
+                      int32_t n,
+                      double* out_scores) {
+    if (!query || qdim <= 0 || !node_ptrs || !node_dims || n <= 0 || !out_scores) return false;
+    if (!eg_metal_ensure_init()) return false;
+
+    @autoreleasepool {
+        const size_t dim = (size_t)qdim;
+        const size_t nu  = (size_t)n;
+
+        /* Gather into a packed row-major matrix — EngramNode.emb is one
+         * malloc per node, not a contiguous array, so this copy is
+         * unavoidable regardless of backend. Rows whose real dim doesn't
+         * match qdim are zero-filled (harmless: the kernel sentinels them
+         * via node_dims before ever reading the row). */
+        float* matrix = (float*)calloc(nu * dim, sizeof(float));
+        int32_t* dims_i32 = (int32_t*)malloc(nu * sizeof(int32_t));
+        if (!matrix || !dims_i32) { free(matrix); free(dims_i32); return false; }
+
+        for (size_t i = 0; i < nu; i++) {
+            dims_i32[i] = node_dims[i];
+            if (node_ptrs[i] && node_dims[i] == qdim) {
+                memcpy(matrix + i * dim, node_ptrs[i], dim * sizeof(float));
+            }
+            /* else: leave zero-filled; node_dims[i] != qdim (or missing)
+             * makes the kernel sentinel it to -2.0 without reading the row. */
+        }
+
+        id<MTLBuffer> bufQuery = [g_device newBufferWithBytes:query
+                                                        length:dim * sizeof(float)
+                                                       options:MTLResourceStorageModeShared];
+        id<MTLBuffer> bufMatrix = [g_device newBufferWithBytes:matrix
+                                                         length:nu * dim * sizeof(float)
+                                                        options:MTLResourceStorageModeShared];
+        id<MTLBuffer> bufDims = [g_device newBufferWithBytes:dims_i32
+                                                       length:nu * sizeof(int32_t)
+                                                      options:MTLResourceStorageModeShared];
+        EgCosineParamsC params = { (uint32_t)nu, (uint32_t)dim };
+        id<MTLBuffer> bufParams = [g_device newBufferWithBytes:&params
+                                                         length:sizeof(params)
+                                                        options:MTLResourceStorageModeShared];
+        id<MTLBuffer> bufOut = [g_device newBufferWithLength:nu * sizeof(float)
+                                                       options:MTLResourceStorageModeShared];
+
+        free(matrix); free(dims_i32);
+
+        if (!bufQuery || !bufMatrix || !bufDims || !bufParams || !bufOut) return false;
+
+        id<MTLCommandBuffer> cmd = [g_queue commandBuffer];
+        if (!cmd) return false;
+        id<MTLComputeCommandEncoder> enc = [cmd computeCommandEncoder];
+        if (!enc) return false;
+
+        [enc setComputePipelineState:g_pipeline];
+        [enc setBuffer:bufQuery  offset:0 atIndex:0];
+        [enc setBuffer:bufMatrix offset:0 atIndex:1];
+        [enc setBuffer:bufDims   offset:0 atIndex:2];
+        [enc setBuffer:bufParams offset:0 atIndex:3];
+        [enc setBuffer:bufOut    offset:0 atIndex:4];
+
+        NSUInteger tgSize = g_pipeline.maxTotalThreadsPerThreadgroup;
+        if (tgSize > 256) tgSize = 256;
+        if (tgSize < 1) tgSize = 1;
+        MTLSize gridSize = MTLSizeMake(nu, 1, 1);
+        MTLSize threadgroupSize = MTLSizeMake(tgSize, 1, 1);
+        [enc dispatchThreads:gridSize threadsPerThreadgroup:threadgroupSize];
+        [enc endEncoding];
+
+        [cmd commit];
+        [cmd waitUntilCompleted];
+
+        if (cmd.status != MTLCommandBufferStatusCompleted) return false;
+
+        const float* results = (const float*)bufOut.contents;
+        if (!results) return false;
+        for (size_t i = 0; i < nu; i++) out_scores[i] = (double)results[i];
+        return true;
+    }
+}
+
+static bool mh_batch_multi(const float* queries, int32_t qdim, int32_t nq,
+                            const float* const* node_ptrs,
+                            const int32_t* node_dims,
+                            int32_t n,
+                            double* out_scores) {
+    if (!queries || qdim <= 0 || nq <= 0 || !node_ptrs || !node_dims || n <= 0 || !out_scores) return false;
+    if (!eg_metal_ensure_init()) return false;
+
+    @autoreleasepool {
+        const size_t dim = (size_t)qdim;
+        const size_t nu  = (size_t)n;
+        const size_t nqu = (size_t)nq;
+
+        float* matrix = (float*)calloc(nu * dim, sizeof(float));
+        int32_t* dims_i32 = (int32_t*)malloc(nu * sizeof(int32_t));
+        if (!matrix || !dims_i32) { free(matrix); free(dims_i32); return false; }
+
+        for (size_t i = 0; i < nu; i++) {
+            dims_i32[i] = node_dims[i];
+            if (node_ptrs[i] && node_dims[i] == qdim) {
+                memcpy(matrix + i * dim, node_ptrs[i], dim * sizeof(float));
+            }
+        }
+
+        /* This is the ONE upload of node_matrix for the whole nq-query batch —
+         * the fix for the measured re-upload-per-query slowdown. */
+        id<MTLBuffer> bufMatrix = [g_device newBufferWithBytes:matrix
+                                                         length:nu * dim * sizeof(float)
+                                                        options:MTLResourceStorageModeShared];
+        id<MTLBuffer> bufDims = [g_device newBufferWithBytes:dims_i32
+                                                       length:nu * sizeof(int32_t)
+                                                      options:MTLResourceStorageModeShared];
+        id<MTLBuffer> bufQueries = [g_device newBufferWithBytes:queries
+                                                          length:nqu * dim * sizeof(float)
+                                                         options:MTLResourceStorageModeShared];
+        EgCosineMultiParamsC params = { (uint32_t)nu, (uint32_t)dim, (uint32_t)nqu };
+        id<MTLBuffer> bufParams = [g_device newBufferWithBytes:&params
+                                                         length:sizeof(params)
+                                                        options:MTLResourceStorageModeShared];
+        id<MTLBuffer> bufOut = [g_device newBufferWithLength:nqu * nu * sizeof(float)
+                                                       options:MTLResourceStorageModeShared];
+
+        free(matrix); free(dims_i32);
+
+        if (!bufMatrix || !bufDims || !bufQueries || !bufParams || !bufOut) return false;
+
+        id<MTLCommandBuffer> cmd = [g_queue commandBuffer];
+        if (!cmd) return false;
+        id<MTLComputeCommandEncoder> enc = [cmd computeCommandEncoder];
+        if (!enc) return false;
+
+        [enc setComputePipelineState:g_pipeline_multi];
+        [enc setBuffer:bufQueries offset:0 atIndex:0];
+        [enc setBuffer:bufMatrix  offset:0 atIndex:1];
+        [enc setBuffer:bufDims    offset:0 atIndex:2];
+        [enc setBuffer:bufParams  offset:0 atIndex:3];
+        [enc setBuffer:bufOut     offset:0 atIndex:4];
+
+        /* 2D dispatch: x over nodes, y over queries. Threadgroup width picked
+         * from the pipeline's own limit, height fixed at 1 — nq is typically
+         * small (tens to low hundreds) relative to n (thousands+), so tiling
+         * the wide axis (n) is what matters for occupancy. */
+        NSUInteger tgWidth = g_pipeline_multi.maxTotalThreadsPerThreadgroup;
+        if (tgWidth > 256) tgWidth = 256;
+        if (tgWidth < 1) tgWidth = 1;
+        MTLSize gridSize = MTLSizeMake(nu, nqu, 1);
+        MTLSize threadgroupSize = MTLSizeMake(tgWidth, 1, 1);
+        [enc dispatchThreads:gridSize threadsPerThreadgroup:threadgroupSize];
+        [enc endEncoding];
+
+        [cmd commit];
+        [cmd waitUntilCompleted];
+
+        if (cmd.status != MTLCommandBufferStatusCompleted) return false;
+
+        const float* results = (const float*)bufOut.contents;
+        if (!results) return false;
+        for (size_t i = 0; i < nqu * nu; i++) out_scores[i] = (double)results[i];
+        return true;
+    }
+}
+
+static const EgCosineBatchStrategy g_metal_hand_strategy = {
+    .name = "metal-hand",
+    .available = mh_available,
+    .batch = mh_batch,
+    .batch_multi = mh_batch_multi,
+};
+
+const EgCosineBatchStrategy* eg_cosine_batch_strategy_metal_hand(void) {
+    return &g_metal_hand_strategy;
+}

lang/runtime/vindex_bench.c

@@ -7,11 +7,29 @@
  *
  * Read-only: never opens a socket, never writes the store. Safe on an nsbx clone.
  *
- * Build:  cc -O2 -std=c11 vindex_bench.c engram_vindex.c -lm -o vindex_bench
+ * Also runs the brute-force oracle a second (and third) way, through the
+ * batch-cosine Strategies behind eg_cosine_batch_strategy.h — the ggml
+ * strategy and the hand-rolled-Metal strategy (Apple/Metal only; see
+ * eg_cosine_batch.h/eg_cosine_batch_strategy.h) — and reports each one's
+ * latency + a correctness check against the CPU oracle side-by-side with the
+ * existing CPU-vs-HNSW numbers. This harness deliberately reaches past the
+ * single-selection Factory (eg_cosine_batch.c) to instantiate every
+ * compiled-in strategy directly, so it can compare all of them against the
+ * SAME dataset in one run — that is the harness's whole job; a real call
+ * site (el_runtime.c) never does this, it only ever calls the plain
+ * eg_cosine_batch()/eg_cosine_batch_multi() adapter functions.
+ * EL_METAL_COSINE=0 forces CPU-only (skips every strategy comparison).
+ *
+ * Build (macOS, ggml + hand-rolled Metal): see build_vindex_bench.sh.
+ * Build (Linux / no Metal): omit every eg_cosine_batch_strategy_*.{c,m} file
+ *   except eg_cosine_batch_strategy_cpu.c — this file never references
+ *   ggml/Metal directly except through the plain-C strategy header, guarded
+ *   by the same EG_HAVE_STRATEGY_* build macros the Factory itself uses.
  * Usage:  vindex_bench store <neuron.egm> <dim> [nqueries] [k] [ef_csv]
  *         vindex_bench synth <N> [dim] [clusters] [nqueries] [k] [ef_csv]
  */
 #include "engram_vindex.h"
+#include "eg_cosine_batch_strategy.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -74,6 +92,106 @@ static double recall_at_k(const int* gt, const uint64_t* ann, int nann, int k){
     return (double)hit / (double)k;
 }
 
+/* EL_METAL_COSINE: 0/off/false disables EVERY strategy comparison outright
+ * (falls back to brute_topk() only), matching el_runtime.c's own gate for
+ * the same env var (back-compat name kept from PR #114; it now gates all
+ * GPU-backed strategies, not just the hand-rolled Metal one). Unset or any
+ * other value = try every compiled-in strategy, report each that's
+ * available, skip (without failing the run) any that isn't. */
+static bool g_strategy_env_checked = false;
+static bool g_strategy_disabled_by_env = false;
+static void eg_strategy_check_env_once(void){
+    if (g_strategy_env_checked) return;
+    g_strategy_env_checked = true;
+    const char* v = getenv("EL_METAL_COSINE");
+    if (v && (v[0]=='0' || v[0]=='n' || v[0]=='N' || v[0]=='f' || v[0]=='F'))
+        g_strategy_disabled_by_env = true;
+}
+
+/* Batched sibling of brute_topk, generalized over ANY EgCosineBatchStrategy:
+ * computes top-k for ALL nq queries in ONE strategy->batch_multi() call,
+ * uploading/preparing the node population exactly once instead of once per
+ * query. out_ids/out_d are nq*k, row-major (query i's results at
+ * out_ids+i*k / out_d+i*k). Returns false (nothing written) on any
+ * failure/unavailability; caller treats that as "skip this strategy in the
+ * report", never as a hard error. */
+static bool batch_topk_strategy(const EgCosineBatchStrategy* strat,
+                                const float* data, int n, int dim,
+                                const float* queries, int nq,
+                                int k, int* out_ids, float* out_d){
+    if (!strat || !strat->available()) return false;
+
+    const float** row_ptrs = malloc((size_t)n * sizeof(float*));
+    int32_t* dims = malloc((size_t)n * sizeof(int32_t));
+    double* scores = malloc((size_t)nq * (size_t)n * sizeof(double));
+    if (!row_ptrs || !dims || !scores) { free(row_ptrs); free(dims); free(scores); return false; }
+
+    for (int i = 0; i < n; i++) { row_ptrs[i] = data + (size_t)i * dim; dims[i] = dim; }
+
+    bool ok = strat->batch_multi(queries, dim, nq, row_ptrs, dims, n, scores);
+    free(row_ptrs); free(dims);
+    if (!ok) { free(scores); return false; }
+
+    for (int qi = 0; qi < nq; qi++) {
+        int* ids = out_ids + (size_t)qi * k;
+        float* ds = out_d + (size_t)qi * k;
+        const double* srow = scores + (size_t)qi * n;
+        for (int i = 0; i < k; i++) { ids[i] = -1; ds[i] = 3.0f; }
+        for (int i = 0; i < n; i++) {
+            float d = 1.0f - (float)srow[i];   /* same distance convention as brute_topk */
+            if (d >= ds[k-1]) continue;
+            int p = k - 1;
+            while (p > 0 && ds[p-1] > d) { ds[p] = ds[p-1]; ids[p] = ids[p-1]; p--; }
+            ds[p] = d; ids[p] = i;
+        }
+    }
+    free(scores);
+    return true;
+}
+
+/* Runs batch_topk_strategy for one named strategy over ALL nq queries, diffs
+ * against the CPU ground truth (gt/gd, both nq*k), and prints a report line
+ * in the same shape PR #114 established for BRUTE-METAL — id-recall over
+ * every query plus the actual max/mean same-rank distance delta across
+ * every (query,rank) pair that was compared, never fabricated or assumed. */
+static void report_strategy_vs_oracle(const char* label, const EgCosineBatchStrategy* strat,
+                                      const float* data, int n, int dim,
+                                      const float* qv, int nq, int k,
+                                      const int* gt, const float* gd, double brute_ms){
+    if (g_strategy_disabled_by_env) { printf("%-13s: disabled via EL_METAL_COSINE\n", label); return; }
+    if (!strat || !strat->available()) { printf("%-13s: not available on this build/host — skipped\n", label); return; }
+
+    int* gtm = malloc((size_t)nq*k*sizeof(int));
+    float* gdm = malloc((size_t)nq*k*sizeof(float));
+    double tm0 = now_s();
+    bool ok = batch_topk_strategy(strat, data, n, dim, qv, nq, k, gtm, gdm);
+    double strat_ms = (now_s()-tm0)*1000.0/nq;
+    if (ok) {
+        double rec_sum = 0; double max_ddiff = 0; double sum_ddiff = 0; int compared = 0;
+        for (int i=0;i<nq;i++) {
+            const int* ids_gt = gt+(size_t)i*k;
+            const float* d_gt = gd+(size_t)i*k;
+            const int* ids_m = gtm+(size_t)i*k;
+            const float* d_m = gdm+(size_t)i*k;
+            uint64_t idset[512]; int m = (k<512)?k:512;
+            for (int j=0;j<m;j++) idset[j] = (uint64_t)ids_m[j];
+            rec_sum += recall_at_k(ids_gt, idset, m, k);
+            for (int j=0;j<k;j++) {
+                if (ids_gt[j] == ids_m[j]) {
+                    double diff = fabs((double)d_gt[j]-(double)d_m[j]);
+                    if (diff>max_ddiff) max_ddiff=diff;
+                    sum_ddiff += diff; compared++;
+                }
+            }
+        }
+        printf("%-13s: %8.3f ms/query   (%.1fx vs CPU brute; id-recall %.4f vs CPU oracle over %d queries; same-rank |Δdist|: max %.2e, mean %.2e over %d compared)\n",
+               label, strat_ms, brute_ms/strat_ms, rec_sum/nq, nq, max_ddiff, compared?sum_ddiff/compared:0.0, compared);
+    } else {
+        printf("%-13s: batch call failed mid-run — skipped\n", label);
+    }
+    free(gtm); free(gdm);
+}
+
 /* Parse "64,128,256" into an int array; returns count. */
 static int parse_csv(const char* s, int* out, int maxo){
     int n=0; if(!s||!*s) return 0;
@@ -142,13 +260,37 @@ static void run_bench(const char* label, float* data, int n, int dim,
         l2norm(dst, dim);
     }
 
-    /* ground truth: brute-force top-k for every query (also the oracle latency). */
+    /* ground truth: brute-force top-k for every query (also the oracle latency).
+     * gd is nq*k (one real slot per query, not a shared scratch buffer) so the
+     * strategy comparisons below can diff against every query's actual
+     * distances, not just whichever query happened to run last. */
     int* gt = malloc((size_t)nq*k*sizeof(int));
-    float* gd = malloc((size_t)k*sizeof(float));
+    float* gd = malloc((size_t)nq*k*sizeof(float));
     double tb0 = now_s();
-    for (int i=0;i<nq;i++) brute_topk(data, n, dim, qv+(size_t)i*dim, k, gt+(size_t)i*k, gd);
+    for (int i=0;i<nq;i++) brute_topk(data, n, dim, qv+(size_t)i*dim, k, gt+(size_t)i*k, gd+(size_t)i*k);
     double brute_ms = (now_s()-tb0)*1000.0/nq;
-    printf("BRUTE-FORCE  : %8.3f ms/query   (oracle; O(N*D))\n", brute_ms);
+    printf("BRUTE-FORCE  : %8.3f ms/query   (oracle; O(N*D), CPU)\n", brute_ms);
+
+    /* GPU-backed oracles: SAME nq queries, SAME top-k contract, via each
+     * compiled-in Strategy's batch_multi() (uploads/prepares the node
+     * population once, not once per query). Run only for strategies that
+     * are actually available (checked internally) — never fabricated, never
+     * assumed. Verified against the CPU ground truth computed above:
+     * id-recall across ALL nq queries, plus the actual max/mean distance
+     * delta across every (query,rank) pair that was compared. */
+    eg_strategy_check_env_once();
+#ifdef EG_HAVE_STRATEGY_GGML
+    report_strategy_vs_oracle("BRUTE-GGML", eg_cosine_batch_strategy_ggml(),
+                              data, n, dim, qv, nq, k, gt, gd, brute_ms);
+#else
+    printf("BRUTE-GGML   : strategy not compiled into this build\n");
+#endif
+#ifdef EG_HAVE_STRATEGY_METAL_HAND
+    report_strategy_vs_oracle("BRUTE-METAL", eg_cosine_batch_strategy_metal_hand(),
+                              data, n, dim, qv, nq, k, gt, gd, brute_ms);
+#else
+    printf("BRUTE-METAL  : strategy not compiled into this build\n");
+#endif
 
     /* HNSW at each ef. */
     uint64_t* aid = malloc((size_t)k*sizeof(uint64_t));