/* 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;
}