diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4da85c41..4c71517d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -205,6 +205,11 @@ llama_build_and_test(test-regex-partial.cpp) llama_build_and_test(test-iq4-ks-kt-decode.cpp) target_include_directories(test-iq4-ks-kt-decode PRIVATE ${PROJECT_SOURCE_DIR}/ggml/src) +# delta-net chunked-recurrence correctness harness (qwen4exp prefill work). +# No model needed: compares sequential ref vs chunk-64 ref vs ggml CPU op. +llama_build(test-delta-chunk.cpp) +llama_test(test-delta-chunk) + # llama_target_and_test(test-opt.cpp) # SLOW llama_target_and_test(test-model-load-cancel.cpp LABEL "model") diff --git a/tests/test-delta-chunk.cpp b/tests/test-delta-chunk.cpp new file mode 100644 index 00000000..0cc1f5bf --- /dev/null +++ b/tests/test-delta-chunk.cpp @@ -0,0 +1,382 @@ +// test-delta-chunk: correctness harness for chunked delta-net (qwen4exp prefill). +// +// Compares three implementations of the fused Gated Delta Rule recurrence +// (see ggml_compute_forward_delta_net_f32 in ggml/src/ggml.c and +// delta_net_recurrent_f32 in ggml/src/ggml-cuda/delta-net.cu): +// (A) plain-C++ sequential reference (exact port, same op order), +// (B) chunked reference (chunk size 64, boundary-state carry) — scaffolding +// for the future parallel-scan CUDA kernel; must match (A) bitwise-ish, +// (C) the repo's ggml op on the CPU backend (exercises the IQK AVX2 path +// on Zen3) — must match (A) within FP reassociation tolerance. +// +// On a CUDA build the same binary gains --cuda mode (TODO): run (C) through +// the CUDA backend and compare vs (A), gating the chunked CUDA kernel. +// +// Usage: test-delta-chunk [--quick] (quick trims the token sweep for CI) + +#include "ggml.h" + +#include +#include +#include +#include +#include +#include + +// ---------------------------------------------------------------- RNG --- +static uint64_t rng_state = 0x123456789abcdefULL; +static float rng_uniform(float lo, float hi) { + rng_state ^= rng_state << 13; + rng_state ^= rng_state >> 7; + rng_state ^= rng_state << 17; + const double u = (double)(rng_state >> 11) * (1.0 / 9007199254740992.0); + return (float)(lo + u * (hi - lo)); +} +static float rng_normal(void) { + const float u1 = rng_uniform(1e-6f, 1.0f); + const float u2 = rng_uniform(0.0f, 1.0f); + return sqrtf(-2.0f * logf(u1)) * cosf(2.0f * (float)M_PI * u2); +} + +// ---------------------------------------------------------------- case --- +struct Case { + int hd; // head dim (64 / 128) + int nt; // n_tokens + int hk; // n_k heads + int hv; // n_v heads + int repeat; // repeat_type 0 (divide) / 1 (mod) + int nseq; // n_seqs + bool saved; // allocate saved_steps + bool hot_state; // large init state (exercises +-1e6 clamp) + bool hot_gate; // large gate values (exercises exp(min(g,50)) cap) + bool loose; // explosive dynamics: relative-tolerance ggml check +}; + +struct Tensors { + std::vector q, k, v, g, beta, state; +}; + +// q,k: [hd,nt,hk,nseq] v: [hd,nt,hv,nseq] state: [hd,hd*hv,nseq] +// g, beta: canonical [nseq,nt,hv] order, idx=(b*nt+t)*hv+h — this matches the +// real graph (build_fused_delta_net permutes scalar gate [hv,t,s] and beta +// [hv,1,t,s] so element (t,h) sits at t*hv+h). The ggml tensor builder below +// transposes into ggml's own memory order; the reference uses canonical idx. +static void fill_tensors(const Case & c, Tensors & t) { + const size_t nq = (size_t)c.hd * c.nt * c.hk * c.nseq; + const size_t nv = (size_t)c.hd * c.nt * c.hv * c.nseq; + const size_t ng = (size_t)c.nt * c.hv * c.nseq; + const size_t ns = (size_t)c.hd * c.hd * c.hv * c.nseq; + t.q.resize(nq); t.k.resize(nq); t.v.resize(nv); + t.g.resize(ng); t.beta.resize(ng); t.state.resize(ns); + for (auto & x : t.q) x = rng_normal(); + for (auto & x : t.k) x = rng_normal(); + // The real graph L2-normalizes q/k before the op (llama-delta-net.cpp), + // and the IQK/CUDA kernels assume pre-normalized inputs — mirror that here + // so the harness tests the kernels, not the (dead in practice) unnorm path. + for (int b = 0; b < c.nseq; ++b) { + for (int h = 0; h < c.hk; ++h) { + for (int tt = 0; tt < c.nt; ++tt) { + float * qq = t.q.data() + ((size_t)b * c.hk + h) * c.hd * c.nt + (size_t)tt * c.hd; + float * kk = t.k.data() + ((size_t)b * c.hk + h) * c.hd * c.nt + (size_t)tt * c.hd; + float qn = 0.0f, kn = 0.0f; + for (int i = 0; i < c.hd; ++i) { qn += qq[i] * qq[i]; kn += kk[i] * kk[i]; } + qn = 1.0f / sqrtf(qn + 1e-12f); + kn = 1.0f / sqrtf(kn + 1e-12f); + for (int i = 0; i < c.hd; ++i) { qq[i] *= qn; kk[i] *= kn; } + } + } + } + for (auto & x : t.v) x = rng_normal() * 0.5f; + // Stable dynamics like the real model (forgetting gates: decay<1 almost + // always). Explosive decay>1 amplifies 1e-7 FP ordering diffs into O(1) + // within a few steps on BOTH sides, which proves nothing — torture cases + // below cover the exp-cap/clamp paths explicitly at small nt. + for (auto & x : t.g) x = c.hot_gate ? rng_uniform(-4.0f, 60.0f) : rng_uniform(-3.0f, -0.1f); + for (auto & x : t.beta) x = rng_uniform(-3.0f, 3.0f); + for (auto & x : t.state) x = c.hot_state ? rng_uniform(-2e6f, 2e6f) : rng_normal() * 0.5f; +} + +// ------------------------------------------------------- sequential ref --- +// Exact port of ggml_compute_forward_delta_net_f32 (fused, scalar-gate path). +// Layouts mirror the kernel (NOT C-order): out is [b,t,h,r] with token stride +// hd*hv (graph views it as [S_v,H_v,nt,nseq]); saved steps are [t,b,h,state] +// with step stride hd*hd*hv*nseq; final state is [b,h,state]. +struct Result { + std::vector out; // [hd,nt,hv,nseq] + std::vector state; // final [hd,hd*hv,nseq] + std::vector saved; // [(nt-1)*hd*hd*hv*nseq] or empty +}; + +static inline float sigmoid_f(float x) { return 1.0f / (1.0f + expf(-x)); } + +static void run_token( + const Case & c, const Tensors & t, int b, int h, int tt, + float * state /*[hd*hd]*/, float * out_t /*[hd]*/, + const float scale) { + const int hd = c.hd; + // NB: repeat_type==1 maps head h to h % H_k (i.e. h % (n_heads/gqa_ratio)), + // NOT h % (H_v/H_k) — matches kernel + ggml_compute_forward_delta_net_f32. + const int hk = (c.repeat == 0) ? h / (c.hv / c.hk) : h % c.hk; + const size_t qoff = ((size_t)b * c.hk + hk) * c.hd * c.nt + (size_t)tt * c.hd; + const size_t voff = ((size_t)b * c.hv + h) * c.hd * c.nt + (size_t)tt * c.hd; + const size_t goff = ((size_t)b * c.nt + tt) * c.hv + h; + const float * q = t.q.data() + qoff; + const float * k = t.k.data() + qoff; + const float * v = t.v.data() + voff; + const float beta = sigmoid_f(t.beta[goff]); + const float decay = expf(fminf(t.g[goff], 50.0f)); + + float qn = 0.0f, kn = 0.0f; + for (int i = 0; i < hd; ++i) { qn += q[i] * q[i]; kn += k[i] * k[i]; } + const float qni = 1.0f / sqrtf(qn + 1e-12f); + const float kni = 1.0f / sqrtf(kn + 1e-12f); + + float score = 0.0f; + for (int i = 0; i < hd; ++i) score += (k[i] * kni) * (q[i] * qni * scale); + + std::vector vn(hd); + for (int r = 0; r < hd; ++r) { + float vp = 0.0f, ov = 0.0f; + for (int col = 0; col < hd; ++col) { + const float s = state[r + col * hd]; + vp += s * k[col]; + ov += s * q[col]; + } + vn[r] = v[r] * beta - vp * beta * decay * kni; + out_t[r] = ov * decay * qni * scale + vn[r] * score; + } + for (int col = 0; col < hd; ++col) { + const float kc = k[col] * kni; + for (int r = 0; r < hd; ++r) { + float s = decay * state[r + col * hd] + vn[r] * kc; + state[r + col * hd] = fminf(fmaxf(s, -1e6f), 1e6f); + } + } +} + +static void run_sequential(const Case & c, const Tensors & t, Result & r) { + const int hd = c.hd; + const float scale = 1.0f / sqrtf((float)hd); + const size_t stsz = (size_t)hd * hd; + r.out.assign((size_t)hd * c.nt * c.hv * c.nseq, 0.0f); + r.state.assign(stsz * c.hv * c.nseq, 0.0f); + r.saved.clear(); + if (c.saved && c.nt > 1) r.saved.assign((size_t)(c.nt - 1) * stsz * c.hv * c.nseq, 0.0f); + std::vector st(stsz); + for (int b = 0; b < c.nseq; ++b) { + for (int h = 0; h < c.hv; ++h) { + memcpy(st.data(), t.state.data() + ((size_t)b * c.hv + h) * stsz, stsz * sizeof(float)); + for (int tt = 0; tt < c.nt; ++tt) { + float * out_t = r.out.data() + (((size_t)b * c.nt + tt) * c.hv + h) * hd; + run_token(c, t, b, h, tt, st.data(), out_t, scale); + if (c.saved && tt + 1 < c.nt) { + memcpy(r.saved.data() + (((size_t)tt * c.nseq + b) * c.hv + h) * stsz, + st.data(), stsz * sizeof(float)); + } + } + memcpy(r.state.data() + ((size_t)b * c.hv + h) * stsz, st.data(), stsz * sizeof(float)); + } + } +} + +// ---------------------------------------------------------- chunked ref --- +// Chunk size 64 (QWEN3NEXT_CHUNK_SIZE). Same per-token math, chunk-at-a-time +// scheduling with explicit boundary-state carry — the structure the future +// parallel-scan kernel must reproduce bit-compatibly. +#define DELTA_CHUNK 64 + +static void run_chunked(const Case & c, const Tensors & t, Result & r) { + const int hd = c.hd; + const float scale = 1.0f / sqrtf((float)hd); + const size_t stsz = (size_t)hd * hd; + r.out.assign((size_t)hd * c.nt * c.hv * c.nseq, 0.0f); + r.state.assign(stsz * c.hv * c.nseq, 0.0f); + r.saved.clear(); + if (c.saved && c.nt > 1) r.saved.assign((size_t)(c.nt - 1) * stsz * c.hv * c.nseq, 0.0f); + std::vector carry(stsz), st(stsz); + for (int b = 0; b < c.nseq; ++b) { + for (int h = 0; h < c.hv; ++h) { + memcpy(carry.data(), t.state.data() + ((size_t)b * c.hv + h) * stsz, stsz * sizeof(float)); + for (int c0 = 0; c0 < c.nt; c0 += DELTA_CHUNK) { + const int c1 = c0 + DELTA_CHUNK < c.nt ? c0 + DELTA_CHUNK : c.nt; + memcpy(st.data(), carry.data(), stsz * sizeof(float)); // chunk-in = prev boundary + for (int tt = c0; tt < c1; ++tt) { + float * out_t = r.out.data() + (((size_t)b * c.nt + tt) * c.hv + h) * hd; + run_token(c, t, b, h, tt, st.data(), out_t, scale); + if (c.saved && tt + 1 < c.nt) { + memcpy(r.saved.data() + (((size_t)tt * c.nseq + b) * c.hv + h) * stsz, + st.data(), stsz * sizeof(float)); + } + } + memcpy(carry.data(), st.data(), stsz * sizeof(float)); // chunk-out boundary + } + memcpy(r.state.data() + ((size_t)b * c.hv + h) * stsz, carry.data(), stsz * sizeof(float)); + } + } +} + +// ------------------------------------------------------------- ggml op --- +static bool run_ggml_op(const Case & c, const Tensors & t, Result & r, std::string & err) { + const int hd = c.hd; + const size_t out_size = (size_t)hd * c.nt * c.hv * c.nseq; + const size_t st_size = (size_t)hd * hd * c.hv * c.nseq; + // result is 1-D [out + final state]; plus headroom ggml may need + const size_t mem = (out_size + st_size + (c.saved && c.nt > 1 ? (size_t)(c.nt - 1) * st_size / (c.hv * c.nseq) * c.hv * c.nseq : 0)) * sizeof(float) + + (t.q.size() + t.k.size() + t.v.size() + t.g.size() + t.beta.size() + t.state.size()) * sizeof(float) + + 64 * 1024 * 1024; + struct ggml_init_params ip = { mem, nullptr, false }; + struct ggml_context * ctx = ggml_init(ip); + if (!ctx) { err = "ggml_init failed"; return false; } + + auto mk = [&](const char * name, ggml_type type, int64_t n0, int64_t n1, int64_t n2, int64_t n3, const void * data, size_t nbytes) { + struct ggml_tensor * ten = ggml_new_tensor_4d(ctx, type, n0, n1, n2, n3); + ggml_set_name(ten, name); + memcpy(ten->data, data, nbytes); + return ten; + }; + struct ggml_tensor * q = mk("q", GGML_TYPE_F32, hd, c.nt, c.hk, c.nseq, t.q.data(), t.q.size() * 4); + struct ggml_tensor * k = mk("k", GGML_TYPE_F32, hd, c.nt, c.hk, c.nseq, t.k.data(), t.k.size() * 4); + struct ggml_tensor * v = mk("v", GGML_TYPE_F32, hd, c.nt, c.hv, c.nseq, t.v.data(), t.v.size() * 4); + // g/beta: the IQK fast path (and the real graph's permuted views) address + // element (t,h) at flat t*hv+h, NOT at the contiguous [nt,1,hv] position + // t+h*nt. Keep canonical [b,t,h] data in place and patch strides to match: + // g [nt,1,hv,nseq]: nb = [hv, nt*hv, 1, nt*hv] (floats). + struct ggml_tensor * g = mk("g", GGML_TYPE_F32, c.nt, 1, c.hv, c.nseq, t.g.data(), t.g.size() * 4); + g->nb[0] = (size_t)c.hv * 4; + g->nb[1] = (size_t)c.nt * c.hv * 4; + g->nb[2] = 4; + g->nb[3] = (size_t)c.nt * c.hv * 4; + struct ggml_tensor * beta = mk("beta", GGML_TYPE_F32, 1, c.nt, c.hv, c.nseq, t.beta.data(), t.beta.size() * 4); + beta->nb[0] = 4; + beta->nb[1] = (size_t)c.hv * 4; + beta->nb[2] = 4; + beta->nb[3] = (size_t)c.nt * c.hv * 4; + struct ggml_tensor * st = mk("st", GGML_TYPE_F32, hd, hd * c.hv, 1, c.nseq, t.state.data(), t.state.size() * 4); + struct ggml_tensor * saved = nullptr; + std::vector saved_buf; + if (c.saved && c.nt > 1) { + saved_buf.assign((size_t)(c.nt - 1) * st_size, 0.0f); + saved = ggml_new_tensor_1d(ctx, GGML_TYPE_F32, (int64_t)saved_buf.size()); + ggml_set_name(saved, "saved"); + memcpy(saved->data, saved_buf.data(), saved_buf.size() * 4); + } + struct ggml_tensor * res = ggml_delta_net(ctx, q, k, v, g, beta, st, saved); + res->op_params[0] = c.repeat; + + struct ggml_cgraph * gf = ggml_new_graph_custom(ctx, 16, false); + ggml_build_forward_expand(gf, res); + if (ggml_graph_compute_with_ctx(ctx, gf, 1) != GGML_STATUS_SUCCESS) { + err = "graph compute failed"; + ggml_free(ctx); + return false; + } + r.out.assign((float *)res->data, (float *)res->data + out_size); + r.state.assign((float *)res->data + out_size, (float *)res->data + out_size + st_size); + if (saved) r.saved.assign((float *)saved->data, (float *)saved->data + saved_buf.size()); + else r.saved.clear(); + ggml_free(ctx); + return true; +} + +// ---------------------------------------------------------------- check --- +static float max_abs_diff(const std::vector & a, const std::vector & b, size_t * at = nullptr) { + float m = 0.0f; + for (size_t i = 0; i < a.size(); ++i) { + if (!std::isfinite(a[i]) || !std::isfinite(b[i])) { if (at) *at = i; return INFINITY; } + const float d = fabsf(a[i] - b[i]); + if (d > m) { m = d; if (at) *at = i; } + } + return m; +} + +// Loose comparator for explosive-dynamics torture cases: same-sign inf counts +// as equal; otherwise relative tolerance on large values, absolute on small. +// Denominator floor is 1e3 (not 1): hot-state intermediates are ~1e5-1e6 scale +// where a single FMA-vs-mul+add ulp is ~1e-2 absolute, and near-cancellation +// can leave small residuals dominated by that noise — all expected FP noise, +// while genuine mapping bugs still show as O(1e3+) and fail loudly. +static float max_rel_diff(const std::vector & a, const std::vector & b, size_t * at = nullptr) { + float m = 0.0f; + for (size_t i = 0; i < a.size(); ++i) { + const float x = a[i], y = b[i]; + if (!std::isfinite(x) || !std::isfinite(y)) { + if (std::isfinite(x) != std::isfinite(y)) { if (at) *at = i; return INFINITY; } + if ((x > 0) != (y > 0)) { if (at) *at = i; return INFINITY; } + continue; + } + const float denom = fabsf(x) > 1e3f ? fabsf(x) : 1e3f; + const float d = fabsf(x - y) / denom; + if (d > m) { m = d; if (at) *at = i; } + } + return m; +} + +static int failures = 0; +static void check_case(const Case & c, int idx) { + Tensors t; + fill_tensors(c, t); + Result seq, chk, op; + run_sequential(c, t, seq); + run_chunked(c, t, chk); + const float d_out_cc = max_abs_diff(seq.out, chk.out); + const float d_st_cc = max_abs_diff(seq.state, chk.state); + float d_sv_cc = 0.0f; + if (c.saved) d_sv_cc = max_abs_diff(seq.saved, chk.saved); + std::string err; + if (!run_ggml_op(c, t, op, err)) { + printf("case %2d hd=%d nt=%d hk=%d hv=%d rep=%d nseq=%d saved=%d hot=%d/%d GGML-OP-FAIL: %s\n", + idx, c.hd, c.nt, c.hk, c.hv, c.repeat, c.nseq, c.saved, c.hot_state, c.hot_gate, err.c_str()); + failures++; + return; + } + const float d_out_op = c.loose ? max_rel_diff(seq.out, op.out) : max_abs_diff(seq.out, op.out); + const float d_st_op = c.loose ? max_rel_diff(seq.state, op.state) : max_abs_diff(seq.state, op.state); + float d_sv_op = 0.0f; + size_t at_sv = 0; + if (c.saved) d_sv_op = c.loose ? max_rel_diff(seq.saved, op.saved, &at_sv) : max_abs_diff(seq.saved, op.saved); + + // chunked must be exact (identical op order); ggml op allows FP reassociation + // (loose torture cases use relative tolerance for explosive dynamics) + const bool ok_cc = d_out_cc <= 1e-6f && d_st_cc <= 1e-6f && d_sv_cc <= 1e-6f; + const float tol = c.loose ? 1e-3f : 2e-4f; + const bool ok_op = d_out_op <= tol && d_st_op <= tol && d_sv_op <= tol; + printf("case %2d hd=%d nt=%3d hk=%d hv=%d rep=%d nseq=%d saved=%d hot=%d/%d " + "chunked[out %.2e st %.2e sv %.2e] %s ggml[out %.2e st %.2e sv %.2e] %s\n", + idx, c.hd, c.nt, c.hk, c.hv, c.repeat, c.nseq, c.saved, c.hot_state, c.hot_gate, + d_out_cc, d_st_cc, d_sv_cc, ok_cc ? "OK " : "FAIL", + d_out_op, d_st_op, d_sv_op, ok_op ? "OK " : "FAIL"); + if (c.loose && !ok_op) { + const size_t stsz = (size_t)c.hd * c.hd; + printf(" worst-sv idx %zu (t=%zu b=%zu h=%zu e=%zu): seq=%.6e ggml=%.6e\n", + at_sv, at_sv / (stsz * c.nseq * c.hv), (at_sv / stsz / c.hv) % c.nseq, + (at_sv / stsz) % c.hv, at_sv % stsz, + seq.saved[at_sv], op.saved[at_sv]); + } + if (!ok_cc || !ok_op) failures++; +} + +int main(int argc, char ** argv) { + const bool quick = argc > 1 && !strcmp(argv[1], "--quick"); + std::vector cases; + const int hds[] = {64, 128}; + const int toks_full[] = {1, 2, 7, 8, 9, 63, 64, 65, 100, 128, 200, 256}; + const int toks_quick[] = {1, 8, 64, 65, 200}; + const int * toks = quick ? toks_quick : toks_full; + const int ntoks = quick ? 5 : 12; + int idx = 0; + for (int hd : hds) { + for (int ti = 0; ti < ntoks; ++ti) { + // (hk,hv): gqa 1 and 4, both repeat types; seqs 1..2; saved on/off + cases.push_back({hd, toks[ti], 4, 4, 0, 1, toks[ti] > 1, false, false, false}); + if (ti % 3 == 0) cases.push_back({hd, toks[ti], 2, 8, 0, 1, toks[ti] > 1, false, false, false}); + if (ti % 3 == 1) cases.push_back({hd, toks[ti], 2, 8, 1, 2, toks[ti] > 1, false, false, false}); + } + } + // torture cases: clamp + exp-cap + saved steps (small nt: hot_gate dynamics + // are explosive by construction, so these pin the paths, not precision) + cases.push_back({128, 9, 2, 8, 0, 1, true, true, true, true}); + cases.push_back({64, 65, 4, 4, 1, 2, true, true, false, true}); + for (auto & c : cases) check_case(c, idx++); + printf("%s: %d/%d cases passed\n", failures ? "FAIL" : "PASS", idx - failures, idx); + return failures ? 1 : 0; +}