ggml_tensor * pe = require_tensor(model, kBase + ".embeddings.position_embeddings");

const int64_t stored_seq = kPretrainGrid * kPretrainGrid; // 576

const size_t n_stored = static_cast<size_t>(stored_seq * kHiddenSize);

std::vector<float> stored(n_stored);

if (pe->type == GGML_TYPE_F32) {

ggml_backend_tensor_get(pe, stored.data(), 0, n_stored * sizeof(float));

} else if (pe->type == GGML_TYPE_F16) {

std::vector<ggml_fp16_t> buf(n_stored);

ggml_backend_tensor_get(pe, buf.data(), 0, n_stored * sizeof(ggml_fp16_t));

for (size_t i = 0; i < n_stored; ++i) stored[i] = ggml_fp16_to_fp32(buf[i]);

} else {

throw std::runtime_error("pos_embed: unsupported dtype");

}

// stored is [576, 1024] in row-major (c fastest). Semantic: [24, 24, 1024].

// Tile to [72, 72, 1024]: for each (h, w) in inference grid, look up

// (h % 24, w % 24) in the pretrained grid.

std::vector<float> tiled(static_cast<size_t>(kGridSize * kGridSize * kHiddenSize));

for (int64_t h = 0; h < kGridSize; ++h) {

for (int64_t w = 0; w < kGridSize; ++w) {

const int64_t ph = h % kPretrainGrid;

const int64_t pw = w % kPretrainGrid;

const size_t src_off = static_cast<size_t>(

(ph * kPretrainGrid + pw) * kHiddenSize);

const size_t dst_off = static_cast<size_t>(

(h * kGridSize + w) * kHiddenSize);

std::memcpy(tiled.data() + dst_off, stored.data() + src_off,

kHiddenSize * sizeof(float));

}

}

return tiled;

const std::string p = layer_prefix(layer);

std::vector<std::string> names;

for (const std::string qkvo : {"q_proj", "k_proj", "v_proj", "o_proj"}) {

names.push_back(p + ".attention." + qkvo + ".weight");

names.push_back(p + ".attention." + qkvo + ".bias");

}

names.push_back(p + ".layer_norm1.weight");

names.push_back(p + ".layer_norm1.bias");

names.push_back(p + ".layer_norm2.weight");

names.push_back(p + ".layer_norm2.bias");

for (const std::string fc : {"fc1", "fc2"}) {

names.push_back(p + ".mlp." + fc + ".weight");

names.push_back(p + ".mlp." + fc + ".bias");

}

return names;

int64_t seq_len; // 576 for window 24×24

int64_t head_dim; // 64

std::vector<float> cos_tab; // [seq_len, head_dim]

std::vector<float> sin_tab; // [seq_len, head_dim]

RopeTable t;

t.head_dim = kHeadDim;

t.seq_len = grid_x * grid_y;

// freqs [dim/4 = 16]

const int64_t nfreqs = kHeadDim / 4;

std::vector<float> freqs(static_cast<size_t>(nfreqs));

for (int64_t i = 0; i < nfreqs; ++i) {

// arange(0, dim, 4) → [0, 4, 8, ..., 60], divide by dim=64

freqs[static_cast<size_t>(i)] = 1.0f / std::pow(kRopeTheta,

static_cast<float>(4 * i) / static_cast<float>(kHeadDim));

}

t.cos_tab.assign(static_cast<size_t>(t.seq_len * t.head_dim), 0.0f);

t.sin_tab.assign(static_cast<size_t>(t.seq_len * t.head_dim), 0.0f);

for (int64_t seq = 0; seq < t.seq_len; ++seq) {

const int64_t x = (seq % grid_x);

const int64_t y = (seq / grid_x);

const float xp = static_cast<float>(x) * scale;

const float yp = static_cast<float>(y) * scale;

// Match PyTorch: inv_freq = cat([outer(x_pos, freqs), outer(y_pos, freqs)])

// inv_freq = repeat_interleave(2, dim=-1)

// So dim k in [0, 32): freq index k/2 with x_pos.

// dim k in [32, 64): freq index (k-32)/2 with y_pos.

for (int64_t k = 0; k < kHeadDim / 2; ++k) {

const float angle = xp * freqs[static_cast<size_t>(k / 2)];

t.cos_tab[static_cast<size_t>(seq * kHeadDim + k)] = std::cos(angle);

t.sin_tab[static_cast<size_t>(seq * kHeadDim + k)] = std::sin(angle);

}

for (int64_t k = kHeadDim / 2; k < kHeadDim; ++k) {

const float angle = yp * freqs[static_cast<size_t>((k - kHeadDim / 2) / 2)];

t.cos_tab[static_cast<size_t>(seq * kHeadDim + k)] = std::cos(angle);

t.sin_tab[static_cast<size_t>(seq * kHeadDim + k)] = std::sin(angle);

}

}

return t;

std::vector<float> & qk, // [B, num_heads, seq, head_dim]

int64_t B, int64_t H, int64_t S, int64_t D,

const RopeTable & tab

) {

if (D != tab.head_dim || S != tab.seq_len) {

throw std::runtime_error("apply_rope_2d_cpu: shape mismatch");

}

for (int64_t b = 0; b < B; ++b) {

for (int64_t h = 0; h < H; ++h) {

for (int64_t s = 0; s < S; ++s) {

float * row = qk.data() + ((b * H + h) * S + s) * D;

const float * c = tab.cos_tab.data() + s * D;

const float * si = tab.sin_tab.data() + s * D;

for (int64_t k = 0; k < D; k += 2) {

const float x0 = row[k];

const float x1 = row[k + 1];

row[k] = x0 * c[k] - x1 * si[k];

row[k + 1] = x1 * c[k + 1] + x0 * si[k + 1];

}

}

}

}

const std::vector<float> & spatial, // [C, W*H, 1] flat, iterating (h, w, c) c fastest

int64_t W, int64_t H, int64_t win

) {

const int64_t nWx = W / win;

const int64_t nWy = H / win;

const int64_t num_windows = nWx * nWy;

std::vector<float> out(spatial.size());

for (int64_t wy = 0; wy < nWy; ++wy) {

for (int64_t wx = 0; wx < nWx; ++wx) {

const int64_t window_idx = wy * nWx + wx;

for (int64_t iy = 0; iy < win; ++iy) {

for (int64_t ix = 0; ix < win; ++ix) {

const int64_t src_h = wy * win + iy;

const int64_t src_w = wx * win + ix;

const int64_t src_seq = src_h * W + src_w;

const int64_t dst_seq = iy * win + ix;

// C is fastest — copy C values contiguously

const size_t src_off = static_cast<size_t>(src_seq) * kHiddenSize;

const size_t dst_off = static_cast<size_t>(

window_idx * win * win + dst_seq) * kHiddenSize;

std::memcpy(out.data() + dst_off, spatial.data() + src_off,

kHiddenSize * sizeof(float));

}

}

}

}

(void)num_windows;

return out;

const std::vector<float> & windowed, // [C, win*win, num_windows] flat

int64_t W, int64_t H, int64_t win

) {

const int64_t nWx = W / win;

const int64_t nWy = H / win;

std::vector<float> out(windowed.size());

for (int64_t wy = 0; wy < nWy; ++wy) {

for (int64_t wx = 0; wx < nWx; ++wx) {

const int64_t window_idx = wy * nWx + wx;

for (int64_t iy = 0; iy < win; ++iy) {

for (int64_t ix = 0; ix < win; ++ix) {

const int64_t dst_h = wy * win + iy;

const int64_t dst_w = wx * win + ix;

const int64_t dst_seq = dst_h * W + dst_w;

const int64_t src_seq = iy * win + ix;

const size_t src_off = static_cast<size_t>(

window_idx * win * win + src_seq) * kHiddenSize;

const size_t dst_off = static_cast<size_t>(dst_seq) * kHiddenSize;

std::memcpy(out.data() + dst_off, windowed.data() + src_off,

kHiddenSize * sizeof(float));

}

}

}

}

return out;

const std::vector<float> & x, int64_t Cin, int64_t H, int64_t W, int64_t Cout,

const std::vector<float> & w, // [Cout, Cin, 1, 1] → flat [Cout, Cin]

const std::vector<float> & b // [Cout]

) {

std::vector<float> y(static_cast<size_t>(Cout * H * W));

const int64_t HW = H * W;

for (int64_t oc = 0; oc < Cout; ++oc) {

const float bo = b[static_cast<size_t>(oc)];

for (int64_t p = 0; p < HW; ++p) {

float s = bo;

for (int64_t ic = 0; ic < Cin; ++ic) {

s += w[static_cast<size_t>(oc * Cin + ic)] *

x[static_cast<size_t>(ic * HW + p)];

}

y[static_cast<size_t>(oc * HW + p)] = s;

}

}

return y;

const std::vector<float> & x, int64_t Cin, int64_t H, int64_t W, int64_t Cout,

const std::vector<float> & w, // [Cout, Cin, 3, 3]

const std::vector<float> & b // [Cout]

) {

std::vector<float> y(static_cast<size_t>(Cout * H * W), 0.0f);

const int64_t K = 3, PAD = 1;

for (int64_t oc = 0; oc < Cout; ++oc) {

const float bo = b[static_cast<size_t>(oc)];

for (int64_t h = 0; h < H; ++h) {

for (int64_t iw = 0; iw < W; ++iw) {

float s = bo;

for (int64_t ic = 0; ic < Cin; ++ic) {

for (int64_t kh = 0; kh < K; ++kh) {

const int64_t src_h = h + kh - PAD;

if (src_h < 0 || src_h >= H) continue;

for (int64_t kw = 0; kw < K; ++kw) {

const int64_t src_w = iw + kw - PAD;

if (src_w < 0 || src_w >= W) continue;

const float xv = x[static_cast<size_t>(

ic * (H * W) + src_h * W + src_w)];

const float wv = w[static_cast<size_t>(

oc * (Cin * K * K) + ic * (K * K) + kh * K + kw)];

s += wv * xv;

}

}

}

y[static_cast<size_t>(oc * (H * W) + h * W + iw)] = s;

}

}

}

return y;

const std::vector<float> & x, int64_t Cin, int64_t H, int64_t W, int64_t Cout,

const std::vector<float> & w, // [Cin, Cout, 2, 2]

const std::vector<float> & b // [Cout]

) {

const int64_t OH = H * 2, OW = W * 2;

std::vector<float> y(static_cast<size_t>(Cout * OH * OW));

for (int64_t oc = 0; oc < Cout; ++oc) {

const float bo = b[static_cast<size_t>(oc)];

for (int64_t out_h = 0; out_h < OH; ++out_h) {

const int64_t in_h = out_h / 2;

const int64_t kh = out_h % 2;

for (int64_t out_w = 0; out_w < OW; ++out_w) {

const int64_t in_w = out_w / 2;

const int64_t kw = out_w % 2;

float s = bo;

for (int64_t ic = 0; ic < Cin; ++ic) {

// weight indexed [ic, oc, kh, kw]

const float wv = w[static_cast<size_t>(

ic * (Cout * 4) + oc * 4 + kh * 2 + kw)];

const float xv = x[static_cast<size_t>(

ic * (H * W) + in_h * W + in_w)];

s += wv * xv;

}

y[static_cast<size_t>(oc * (OH * OW) + out_h * OW + out_w)] = s;

}

}

}

return y;

return {

kBase + ".embeddings.patch_embeddings.projection.weight",

kBase + ".embeddings.position_embeddings",

kBase + ".layer_norm.weight",

kBase + ".layer_norm.bias",

};

size_t probed = 0;

for (const auto & name : trunk_top_tensor_names()) {

if (model_.find_tensor(name) == nullptr) {

throw std::runtime_error("vision_encoder: missing trunk-top tensor: " + name);

}

++probed;

}

for (int layer = 0; layer < kNumLayers; ++layer) {

for (const auto & name : per_layer_tensor_names(layer)) {

if (model_.find_tensor(name) == nullptr) {

throw std::runtime_error("vision_encoder: missing per-layer tensor: " + name);

}

++probed;

}

}

return probed;

const std::vector<float> & pixel_values,

const std::vector<int64_t> & pixel_shape

) const {

if (pixel_shape.size() != 3 || pixel_shape[0] != 3 ||

pixel_shape[1] != kImageSize || pixel_shape[2] != kImageSize) {

throw std::runtime_error("pixel_values must be [3, 1008, 1008]");

}

ggml_backend_t backend = model_.backend();

ggml_backend_t cpu_backend = nullptr;

ggml_backend_t backends[2] = { backend, nullptr };

int n_backends = 1;

if (ggml_backend_dev_type(ggml_backend_get_device(backend)) != GGML_BACKEND_DEVICE_TYPE_CPU) {

cpu_backend = model_.cpu_backend();

if (cpu_backend != nullptr) { backends[1] = cpu_backend; n_backends = 2; }

}

const size_t graph_size = 4096;

const size_t ctx_size = ggml_tensor_overhead() * graph_size +

ggml_graph_overhead_custom(graph_size, false);

std::vector<uint8_t> ctx_buf(ctx_size);

ggml_context * ctx = ggml_init({ctx_size, ctx_buf.data(), true});

ggml_backend_sched_t sched = ggml_backend_sched_new(backends, nullptr, n_backends, graph_size, false, true);

ggml_cgraph * gf = ggml_new_graph_custom(ctx, graph_size, false);

ggml_tensor * pv_t = ggml_new_tensor_4d(ctx, GGML_TYPE_F32, kImageSize, kImageSize, 3, 1);

ggml_backend_sched_set_tensor_backend(sched, pv_t, backend);

// No weight permute: PyTorch [OutC, InC, KH, KW] → GGUF ne=[KW, KH, InC, OutC]

// exactly matches ggml_conv_2d's expected [KW, KH, IC, OC] layout.

ggml_tensor * proj_w = ensure_f32(ctx, require_tensor(model_,

kBase + ".embeddings.patch_embeddings.projection.weight"));

ggml_tensor * conv = ggml_conv_2d(ctx, proj_w, pv_t,

kPatchSize, kPatchSize, 0, 0, 1, 1);

// conv ne = [W=72, H=72, C=1024, 1]. Match PyTorch flatten(2).transpose(1,2)

// = [N, H*W, C]: permute so C is fastest, then W, then H → matches

// PyTorch memory layout after the transpose.

// ggml_permute uses `result.ne[axis_i] = a.ne[i]` — so `permute(a, 1, 2, 0, 3)`

// takes conv output ne=[W, H, C, 1] and produces new ne=[C, W, H, 1] with

// C fastest, W middle, H outer. Materialized memory order: (h, w, c) with

// c innermost — exactly matches PyTorch's flatten(2).transpose(1,2) layout.

ggml_tensor * out = ggml_cont(ctx, ggml_permute(ctx, conv, 1, 2, 0, 3));

out = ggml_reshape_3d(ctx, out, kHiddenSize, kNumPatches, 1);

ggml_tensor * cap = ggml_cont(ctx, out);

ggml_build_forward_expand(gf, cap);

ggml_backend_sched_alloc_graph(sched, gf);

ggml_backend_tensor_set(pv_t, pixel_values.data(), 0, pixel_values.size() * sizeof(float));

if (ggml_backend_sched_graph_compute(sched, gf) != GGML_STATUS_SUCCESS) {

ggml_backend_sched_free(sched); ggml_free(ctx);

throw std::runtime_error("patch_embed compute failed");

}

std::vector<float> result(static_cast<size_t>(ggml_nelements(cap)));

ggml_backend_tensor_get(cap, result.data(), 0, result.size() * sizeof(float));

ggml_backend_sched_free(sched);

ggml_free(ctx);

return result;

const std::vector<float> & pixel_values,

const std::vector<int64_t> & pixel_shape

) const {

if (pixel_shape.size() != 3 ||

pixel_shape[0] != 3 ||

pixel_shape[1] != kImageSize ||

pixel_shape[2] != kImageSize) {

throw std::runtime_error("vision_encoder: pixel_values must be [3, 1008, 1008]");

}

ggml_backend_t backend = model_.backend();

ggml_backend_t cpu_backend = nullptr;

ggml_backend_t backends[2] = { backend, nullptr };

int n_backends = 1;

if (ggml_backend_dev_type(ggml_backend_get_device(backend)) != GGML_BACKEND_DEVICE_TYPE_CPU) {

cpu_backend = model_.cpu_backend();

if (cpu_backend != nullptr) { backends[1] = cpu_backend; n_backends = 2; }

}

const size_t graph_size = 8192;

const size_t ctx_size = ggml_tensor_overhead() * graph_size +

ggml_graph_overhead_custom(graph_size, false);

std::vector<uint8_t> ctx_buf(ctx_size);

ggml_context * ctx = ggml_init({ctx_size, ctx_buf.data(), true});

if (!ctx) throw std::runtime_error("ggml_init failed");

ggml_backend_sched_t sched = ggml_backend_sched_new(backends, nullptr, n_backends, graph_size, false, true);

if (!sched) { ggml_free(ctx); throw std::runtime_error("sched_new failed"); }

ggml_cgraph * gf = ggml_new_graph_custom(ctx, graph_size, false);

// pixel_values in ggml column-major: ne = [W, H, C, N] = [1008, 1008, 3, 1]

ggml_tensor * pv_t = ggml_new_tensor_4d(ctx, GGML_TYPE_F32,

kImageSize, kImageSize, 3, 1);

ggml_backend_sched_set_tensor_backend(sched, pv_t, backend);

// Patch embed: Conv2d(3→1024, k=14, s=14, no bias)

// PyTorch weight [Out=1024, In=3, KH=14, KW=14] → GGUF ne [14, 14, 3, 1024]

// = exactly what ggml_conv_2d expects [KW, KH, InC, OutC].

ggml_tensor * proj_w = ensure_f32(ctx, require_tensor(model_,

kBase + ".embeddings.patch_embeddings.projection.weight"));

ggml_tensor * conv = ggml_conv_2d(ctx, proj_w, pv_t,

kPatchSize, kPatchSize, // stride

0, 0, // padding

1, 1); // dilation

// conv now has ne = [72, 72, 1024, 1]

// Reshape spatial [W=72, H=72, C=1024, N=1] to token layout

// [C=1024, seq=5184, 1, 1] — matches Sam3ViTEmbeddings.forward's

// flatten(2).transpose(1,2): [B, C, H, W] → [B, C, H*W] → [B, H*W, C].

ggml_tensor * tokens = ggml_reshape_3d(ctx, conv, kHiddenSize, kNumPatches, 1);

// Now tokens ne = [1024, 5184, 1] which in row-major is [seq, hidden] with

// dim 0 = hidden (fastest). BUT ggml views are column-major-first, so this

// shape means "hidden is fastest, then seq". PyTorch semantically wants

// [seq_len, hidden] where hidden is trailing/fastest. So this IS correct.

// Wait — need to be more careful. conv output has spatial layout

// [W, H, C, 1]. Flatten(2) in PyTorch reduces [B, C, H, W] to

// [B, C, H*W]. That's C outer, H*W inner. Then transpose(1,2) swaps to

// [B, H*W, C]. So end result is [seq, C] with C fastest.

// For ggml (column-major), [seq, C] as PyTorch view corresponds to

// ne=[C, seq] where C is fastest. That matches ggml_reshape_3d(conv, C, seq, 1).

// BUT — conv output layout in ggml: ne=[W, H, C, 1], stride=[1, W, W*H, ...].

// W is fastest. To get PyTorch [B, C, H*W] layout, we need to permute

// dim 0 (W) with dim 2 (C) so C becomes fastest. Then H*W folds.

//

// Actually the sam3cpp vision_trunk.cpp faces the exact same problem

// and handles it. Let me use the same pattern: transpose spatial dims

// together to get channel-fastest layout.

//

// For now: use ggml_cont(ggml_permute(...)) to bring C to fastest,

// then reshape to [C, HW].

// Corrected permute: ggml_permute uses `result.ne[axis_i] = a.ne[i]`.

// permute(conv, 1, 2, 0, 3) with conv ne=[W, H, C, 1] gives result

// ne=[C, W, H, 1] with C fastest, W middle, H outer — matches PyTorch's

// flatten(2).transpose(1, 2) memory layout.

tokens = ggml_cont(ctx, ggml_permute(ctx, conv, 1, 2, 0, 3));

tokens = ggml_reshape_3d(ctx, tokens, kHiddenSize, kNumPatches, 1);

// Add position embeddings — tiled from [1, 576, 1024] (24×24 pretrain

// grid) to [1, 5184, 1024] (72×72 inference grid). Tiling is CPU-side.

const std::vector<float> tiled_pos = load_and_tile_pos_embed(model_);

ggml_tensor * pos_embed_t = ggml_new_tensor_3d(ctx, GGML_TYPE_F32,

kHiddenSize, kNumPatches, 1);

ggml_backend_sched_set_tensor_backend(sched, pos_embed_t, backend);

tokens = ggml_add(ctx, tokens, pos_embed_t);

// Pre-trunk LayerNorm

tokens = layer_norm(ctx, tokens,

require_tensor(model_, kBase + ".layer_norm.weight"),

require_tensor(model_, kBase + ".layer_norm.bias"));

// Reshape to spatial [W, H, C, 1] for windowed attention downstream —

// but for validation, keep as [C, seq, 1] since PyTorch reference stores

// this intermediate in the same format.

ggml_tensor * out = ggml_cont(ctx, tokens);

ggml_build_forward_expand(gf, out);

ggml_backend_sched_alloc_graph(sched, gf);

ggml_backend_tensor_set(pv_t, pixel_values.data(), 0, pixel_values.size() * sizeof(float));

ggml_backend_tensor_set(pos_embed_t, tiled_pos.data(), 0, tiled_pos.size() * sizeof(float));

if (ggml_backend_sched_graph_compute(sched, gf) != GGML_STATUS_SUCCESS) {

ggml_backend_sched_free(sched); ggml_free(ctx);

throw std::runtime_error("vision_encoder prenorm compute failed");

}

std::vector<float> result(static_cast<size_t>(ggml_nelements(out)));

ggml_backend_tensor_get(out, result.data(), 0, result.size() * sizeof(float));

ggml_backend_sched_free(sched);

ggml_free(ctx);

return result;

const std::vector<float> & pixel_values,

const std::vector<int64_t> & pixel_shape

) const {

auto hidden = run_prenorm(pixel_values, pixel_shape);

auto tstart = std::chrono::steady_clock::now();

for (int layer = 0; layer < kNumLayers; ++layer) {

auto tl = std::chrono::steady_clock::now();

hidden = run_layer(layer, hidden);

const auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(

std::chrono::steady_clock::now() - tl).count();

fprintf(stderr, " [vis] layer %d done (%lld ms)\n", layer, (long long)ms); fflush(stderr);

}

return hidden;

int layer_idx, const std::vector<float> & hidden_in

) const {

if (hidden_in.size() != static_cast<size_t>(kNumPatches * kHiddenSize)) {

throw std::runtime_error("run_layer: expected [72*72, 1024] input");

}

const std::string p = layer_prefix(layer_idx);

const bool global = is_global_layer(layer_idx);

auto get_f32 = [&](const std::string & name, size_t n) {

ggml_tensor * t = require_tensor(model_, name);

std::vector<float> v(n);

if (t->type == GGML_TYPE_F32) {

ggml_backend_tensor_get(t, v.data(), 0, v.size() * sizeof(float));

} else if (t->type == GGML_TYPE_F16) {

std::vector<ggml_fp16_t> buf(n);

ggml_backend_tensor_get(t, buf.data(), 0, buf.size() * sizeof(ggml_fp16_t));

for (size_t i = 0; i < n; ++i) v[i] = ggml_fp16_to_fp32(buf[i]);

} else throw std::runtime_error("run_layer: unsupported dtype");

return v;

};

// Same pattern as the LM: matmul weights are F16 in the GGUF, so read

// them AS F16 straight through — ggml_mul_mat's F16×F32 SIMD kernel

// does the dequant per-lane and halves the weight-streaming bandwidth

// vs get_f32.

auto get_f16_raw = [&](const std::string & name, size_t n) {

ggml_tensor * t = require_tensor(model_, name);

if (t->type != GGML_TYPE_F16) {

throw std::runtime_error("get_f16_raw: expected F16 for " + name);

}

std::vector<ggml_fp16_t> v(n);

ggml_backend_tensor_get(t, v.data(), 0, v.size() * sizeof(ggml_fp16_t));

return v;

};

const auto ln1_w = get_f32(p + ".layer_norm1.weight", kHiddenSize);

const auto ln1_b = get_f32(p + ".layer_norm1.bias", kHiddenSize);

const auto ln2_w = get_f32(p + ".layer_norm2.weight", kHiddenSize);

const auto ln2_b = get_f32(p + ".layer_norm2.bias", kHiddenSize);

const std::vector<float> residual1 = hidden_in;

auto cpu_ln = [&](const std::vector<float> & x, int64_t N, int64_t D,

const std::vector<float> & w, const std::vector<float> & b) {

std::vector<float> out(x.size());

for (int64_t i = 0; i < N; ++i) {

double mean = 0.0, var = 0.0;

for (int64_t d = 0; d < D; ++d) mean += x[static_cast<size_t>(i * D + d)];

mean /= D;

for (int64_t d = 0; d < D; ++d) {

const double diff = x[static_cast<size_t>(i * D + d)] - mean;

var += diff * diff;

}

var /= D;

const double inv_std = 1.0 / std::sqrt(var + kLayerNormEps);

for (int64_t d = 0; d < D; ++d) {

const double normed = (x[static_cast<size_t>(i * D + d)] - mean) * inv_std;

out[static_cast<size_t>(i * D + d)] = static_cast<float>(

normed * w[static_cast<size_t>(d)] + b[static_cast<size_t>(d)]);

}

}

return out;

};

const std::vector<float> ln1_out = cpu_ln(hidden_in, kNumPatches, kHiddenSize, ln1_w, ln1_b);

const int64_t win = global ? kGridSize : kWindowSize;

const int64_t nW = global ? 1 : (kGridSize / kWindowSize);

const int64_t num_windows = nW * nW;

const int64_t seq_per_win = win * win;

std::vector<float> windowed = global

? ln1_out

: cpu_window_partition(ln1_out, kGridSize, kGridSize, win);

const auto qw = get_f16_raw(p + ".attention.q_proj.weight", kHiddenSize * kHiddenSize);

const auto qb = get_f32(p + ".attention.q_proj.bias", kHiddenSize);

const auto kw = get_f16_raw(p + ".attention.k_proj.weight", kHiddenSize * kHiddenSize);

const auto kb = get_f32(p + ".attention.k_proj.bias", kHiddenSize);

const auto vw = get_f16_raw(p + ".attention.v_proj.weight", kHiddenSize * kHiddenSize);

const auto vb = get_f32(p + ".attention.v_proj.bias", kHiddenSize);

const auto ow = get_f16_raw(p + ".attention.o_proj.weight", kHiddenSize * kHiddenSize);

const auto ob = get_f32(p + ".attention.o_proj.bias", kHiddenSize);

using detail::cpu_linear;

const int64_t total_tokens = num_windows * seq_per_win;

auto Q = cpu_linear(windowed, total_tokens, kHiddenSize, kHiddenSize, qw, nullptr);

#pragma omp parallel for schedule(static)

for (int64_t n = 0; n < total_tokens; ++n)

for (int64_t d = 0; d < kHiddenSize; ++d)

Q[n * kHiddenSize + d] += qb[d];

auto K = cpu_linear(windowed, total_tokens, kHiddenSize, kHiddenSize, kw, nullptr);

#pragma omp parallel for schedule(static)

for (int64_t n = 0; n < total_tokens; ++n)

for (int64_t d = 0; d < kHiddenSize; ++d)

K[n * kHiddenSize + d] += kb[d];

auto V = cpu_linear(windowed, total_tokens, kHiddenSize, kHiddenSize, vw, nullptr);

#pragma omp parallel for schedule(static)

for (int64_t n = 0; n < total_tokens; ++n)

for (int64_t d = 0; d < kHiddenSize; ++d)

V[n * kHiddenSize + d] += vb[d];

// Transpose to [num_windows, num_heads, seq_per_win, head_dim] for RoPE + attention

auto to_heads = [&](const std::vector<float> & flat) {

std::vector<float> out(flat.size());

for (int64_t b = 0; b < num_windows; ++b) {

for (int64_t h = 0; h < kNumHeads; ++h) {

for (int64_t s = 0; s < seq_per_win; ++s) {

for (int64_t d = 0; d < kHeadDim; ++d) {

const size_t src = ((b * seq_per_win + s) * kNumHeads + h) * kHeadDim + d;

const size_t dst = ((b * kNumHeads + h) * seq_per_win + s) * kHeadDim + d;

out[dst] = flat[src];

}

}

}

}

return out;

};

Q = to_heads(Q); K = to_heads(K); V = to_heads(V);

const float rope_scale = global ? (static_cast<float>(kWindowSize) / kGridSize) : 1.0f;

const RopeTable rope = build_rope_table_2d(win, win, rope_scale);

apply_rope_2d_cpu(Q, num_windows, kNumHeads, seq_per_win, kHeadDim, rope);

apply_rope_2d_cpu(K, num_windows, kNumHeads, seq_per_win, kHeadDim, rope);

// Scaled dot-product attention via ggml_flash_attn_ext.

// Input Q/K/V layout: (num_windows, kNumHeads, seq_per_win, kHeadDim)

// with kHeadDim innermost — matches the (batch, n_head, seq, head_dim)

// memory shape flash_attn expects.

const float attn_scale = 1.0f / std::sqrt(static_cast<float>(kHeadDim));

std::vector<float> attn_hns = detail::flash_attn(

Q, K, V,

/*batch=*/num_windows, /*n_head=*/kNumHeads, /*n_head_kv=*/kNumHeads,

/*seq_q=*/seq_per_win, /*seq_kv=*/seq_per_win, /*head_dim=*/kHeadDim,

attn_scale);

// flash_attn returns (batch, n_head, seq, head_dim). Repack into

// (batch, seq, n_head*head_dim) = attn[total_tokens, kHiddenSize] to

// match what the O-projection expects, mirroring the old code's tail.

std::vector<float> attn(static_cast<size_t>(total_tokens * kHiddenSize));

#pragma omp parallel for schedule(static)

for (int64_t b = 0; b < num_windows; ++b) {

for (int64_t h = 0; h < kNumHeads; ++h) {

for (int64_t s = 0; s < seq_per_win; ++s) {

const size_t src = ((b * kNumHeads + h) * seq_per_win + s) * kHeadDim;

const size_t dst = ((b * seq_per_win + s) * kNumHeads + h) * kHeadDim;

std::memcpy(attn.data() + dst, attn_hns.data() + src,

kHeadDim * sizeof(float));

}

}

}

auto attn_proj = cpu_linear(attn, total_tokens, kHiddenSize, kHiddenSize, ow, nullptr);

for (int64_t n = 0; n < total_tokens; ++n)

for (int64_t d = 0; d < kHiddenSize; ++d)

attn_proj[n * kHiddenSize + d] += ob[d];

std::vector<float> attn_unwin = global

? std::move(attn_proj)

: cpu_window_unpartition(attn_proj, kGridSize, kGridSize, win);

for (size_t i = 0; i < attn_unwin.size(); ++i) attn_unwin[i] += residual1[i];

const std::vector<float> residual2 = attn_unwin;

const std::vector<float> ln2_out = cpu_ln(attn_unwin, kNumPatches, kHiddenSize, ln2_w, ln2_b);

const auto fc1_w = get_f16_raw(p + ".mlp.fc1.weight", kMlpDim * kHiddenSize);

const auto fc1_b = get_f32 (p + ".mlp.fc1.bias", kMlpDim);

const auto fc2_w = get_f16_raw(p + ".mlp.fc2.weight", kHiddenSize * kMlpDim);

const auto fc2_b = get_f32 (p + ".mlp.fc2.bias", kHiddenSize);

auto mlp_mid = cpu_linear(ln2_out, kNumPatches, kHiddenSize, kMlpDim, fc1_w, nullptr);

for (int64_t n = 0; n < kNumPatches; ++n)

for (int64_t d = 0; d < kMlpDim; ++d)

mlp_mid[n * kMlpDim + d] += fc1_b[d];

for (float & x : mlp_mid) {

x = 0.5f * x * (1.0f + std::erf(x / std::sqrt(2.0f)));

}

auto mlp_out = cpu_linear(mlp_mid, kNumPatches, kMlpDim, kHiddenSize, fc2_w, nullptr);

for (int64_t n = 0; n < kNumPatches; ++n)

for (int64_t d = 0; d < kHiddenSize; ++d)

mlp_out[n * kHiddenSize + d] += fc2_b[d];

for (size_t i = 0; i < mlp_out.size(); ++i) mlp_out[i] += residual2[i];

return mlp_out;

const std::vector<float> & trunk_output

) const {

if (trunk_output.size() != static_cast<size_t>(kNumPatches * kHiddenSize)) {

throw std::runtime_error("run_neck: trunk_output must be [72*72, 1024]");

}

// Trunk output is [seq=5184, hidden=1024] with hidden fastest, i.e.

// memory order (h, w, c) with c inner. FPN expects [C, H, W] with

// w fastest. Reshape: for each (c, h, w), target[c*72*72 + h*72 + w]

// = trunk_output[(h*72 + w) * 1024 + c].

std::vector<float> hs(static_cast<size_t>(kHiddenSize * kGridSize * kGridSize));

for (int64_t h = 0; h < kGridSize; ++h) {

for (int64_t w = 0; w < kGridSize; ++w) {

for (int64_t c = 0; c < kHiddenSize; ++c) {

hs[static_cast<size_t>(c * kGridSize * kGridSize + h * kGridSize + w)] =

trunk_output[static_cast<size_t>((h * kGridSize + w) * kHiddenSize + c)];

}

}

}

const std::string base = "backbone.vision_backbone.convs.fpn_layers";

auto get_f32 = [&](const std::string & name, size_t n) {

ggml_tensor * t = require_tensor(model_, name);

std::vector<float> v(n);

if (t->type == GGML_TYPE_F32) {

ggml_backend_tensor_get(t, v.data(), 0, v.size() * sizeof(float));

} else if (t->type == GGML_TYPE_F16) {

std::vector<ggml_fp16_t> buf(n);

ggml_backend_tensor_get(t, buf.data(), 0, buf.size() * sizeof(ggml_fp16_t));

for (size_t i = 0; i < n; ++i) v[i] = ggml_fp16_to_fp32(buf[i]);

} else throw std::runtime_error("run_neck: unsupported dtype");

return v;

};

FpnOutputs out;

// Native ggml FPN neck. Each conv uses ggml_conv_2d (or transpose_2d)

// with pre-computed weights loaded once via get_f32. Bias adds are OMP

// loops after the conv to avoid the ggml_add-broadcast crash at large N.

auto convT_ks = [&](const std::vector<float> & x, int64_t Cin, int64_t H, int64_t W,

int64_t Cout, const std::vector<float> & w,

const std::vector<float> & b, int stride) {

return detail::ggml_conv_transpose2d_ks(x, Cin, H, W, Cout, /*KH=*/2, /*KW=*/2,

GGML_TYPE_F32, w.data(), &b, stride);

};

auto conv1x1 = [&](const std::vector<float> & x, int64_t Cin, int64_t H, int64_t W,

int64_t Cout, const std::vector<float> & w,

const std::vector<float> & b) {

return detail::ggml_conv2d(x, Cin, H, W, Cout, /*KH=*/1, /*KW=*/1,

GGML_TYPE_F32, w.data(), &b, /*stride=*/1, /*pad=*/0);

};

auto conv3x3 = [&](const std::vector<float> & x, int64_t Cin, int64_t H, int64_t W,

int64_t Cout, const std::vector<float> & w,

const std::vector<float> & b) {

return detail::ggml_conv2d(x, Cin, H, W, Cout, /*KH=*/3, /*KW=*/3,

GGML_TYPE_F32, w.data(), &b, /*stride=*/1, /*pad=*/1);

};

auto gelu_inplace = [](std::vector<float> & v) {

#pragma omp parallel for schedule(static)

for (int64_t i = 0; i < static_cast<int64_t>(v.size()); ++i) {

v[i] = 0.5f * v[i] * (1.0f + std::erf(v[i] / std::sqrt(2.0f)));

}

};

// ── FPN layer 0 (scale=4×): ConvT(1024→512) → GELU → ConvT(512→256)

// → proj1 Conv1x1(256→256) → proj2 Conv3x3(256→256) ─

{

const std::string p = base + ".0";

const auto sl0_w = get_f32(p + ".scale_layers.0.weight", 1024 * 512 * 4);

const auto sl0_b = get_f32(p + ".scale_layers.0.bias", 512);

const auto sl2_w = get_f32(p + ".scale_layers.2.weight", 512 * 256 * 4);

const auto sl2_b = get_f32(p + ".scale_layers.2.bias", 256);

const auto p1_w = get_f32(p + ".proj1.weight", 256 * 256);

const auto p1_b = get_f32(p + ".proj1.bias", 256);

const auto p2_w = get_f32(p + ".proj2.weight", 256 * 256 * 9);

const auto p2_b = get_f32(p + ".proj2.bias", 256);

auto t = convT_ks(hs, 1024, kGridSize, kGridSize, 512, sl0_w, sl0_b, 2);

gelu_inplace(t);

t = convT_ks(t, 512, kGridSize * 2, kGridSize * 2, 256, sl2_w, sl2_b, 2);

t = conv1x1(t, 256, kGridSize * 4, kGridSize * 4, 256, p1_w, p1_b);

out.bb0 = conv3x3(t, 256, kGridSize * 4, kGridSize * 4, 256, p2_w, p2_b);

}

// ── FPN layer 1 (scale=2×): ConvT(1024→512) → proj1 (512→256) → proj2 (256→256) ─

{

const std::string p = base + ".1";

const auto sl0_w = get_f32(p + ".scale_layers.0.weight", 1024 * 512 * 4);

const auto sl0_b = get_f32(p + ".scale_layers.0.bias", 512);

const auto p1_w = get_f32(p + ".proj1.weight", 256 * 512);

const auto p1_b = get_f32(p + ".proj1.bias", 256);

const auto p2_w = get_f32(p + ".proj2.weight", 256 * 256 * 9);

const auto p2_b = get_f32(p + ".proj2.bias", 256);

auto t = convT_ks(hs, 1024, kGridSize, kGridSize, 512, sl0_w, sl0_b, 2);

t = conv1x1(t, 512, kGridSize * 2, kGridSize * 2, 256, p1_w, p1_b);

out.bb1 = conv3x3(t, 256, kGridSize * 2, kGridSize * 2, 256, p2_w, p2_b);

}

// ── FPN layer 2 (scale=1×): proj1 (1024→256) → proj2 (256→256) ────────

{

const std::string p = base + ".2";

const auto p1_w = get_f32(p + ".proj1.weight", 256 * 1024);

const auto p1_b = get_f32(p + ".proj1.bias", 256);

const auto p2_w = get_f32(p + ".proj2.weight", 256 * 256 * 9);

const auto p2_b = get_f32(p + ".proj2.bias", 256);

auto t = conv1x1(hs, 1024, kGridSize, kGridSize, 256, p1_w, p1_b);

out.bb2 = conv3x3(t, 256, kGridSize, kGridSize, 256, p2_w, p2_b);

}

return out;