/* engram_vindex.c — HNSW ANN index over f32 embedding vectors (design §9 M8). * * Self-contained: plain C11, stdlib + libm (-lm for sqrtf/logf) only. No * dependency on el_runtime; the store is read via its PERMANENT on-disk format * (design §2.4), decoded read-only here so engram_store.{c,h} stay untouched. * * Algorithm: Malkov & Yashunin, "Efficient and robust approximate nearest * neighbor search using Hierarchical Navigable Small World graphs" (2016). * - multi-layer graph; level ~ Exp(1/ln M), assigned by a per-node seeded PRNG * (deterministic: seed = FIXED_SEED ^ node_ordinal) so a rebuild is bit-for- * bit reproducible regardless of wall-clock or global rand() state. * - greedy descent through upper layers to an entry point, then an ef-bounded * best-first search at each layer (Algorithm 2). * - neighbour selection by the diversity heuristic (Algorithm 4), not plain * k-nearest, with keep-pruned backfill for connectivity. * - bidirectional links; a neighbour whose degree exceeds M (2M on layer 0) is * re-pruned with the same heuristic. * * Metric: vectors are L2-normalised on entry, so cosine similarity == dot * product; distance = 1 - dot (in [0,2], smaller == nearer). Deterministic tie- * breaks are by element index so results are stable across identical builds. */ #include "engram_vindex.h" #include #include #include #include #include #include #include #include /* Deterministic PRNG seed base (fixed constant — never wall-clock/rand). */ #define VINDEX_FIXED_SEED 0x9E3779B97F4A7C15ULL /* ── deterministic PRNG (splitmix64) ──────────────────────────────────────── */ static inline uint64_t splitmix64(uint64_t* s){ uint64_t z = (*s += 0x9E3779B97F4A7C15ULL); z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL; z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL; return z ^ (z >> 31); } /* Uniform double in (0,1]. */ static inline double sm_uniform(uint64_t* s){ /* 53-bit mantissa; +1 keeps it in (0,1] so log() never sees 0. */ return ((double)((splitmix64(s) >> 11) + 1)) * (1.0 / 9007199254740993.0); } /* ── element + index structures ───────────────────────────────────────────── */ typedef struct { int count; int cap; int* ids; /* neighbour element indices */ } NeighList; typedef struct { uint64_t node_id; int level; /* top layer this element appears on (>=0) */ float* vec; /* dim floats, L2-normalised */ NeighList* links; /* level+1 lists; links[l] = neighbours at layer l */ } Elem; struct VIndex { int dim; int M; /* max neighbours per node, upper layers */ int M0; /* == 2*M, layer 0 */ int ef_construction; double mL; /* level normaliser = 1/ln(M) */ Elem* elems; size_t n; size_t cap; int entry; /* entry-point element index, -1 if empty */ int max_level; /* current top layer */ /* scratch: version-stamped visited set (O(1) reset). */ uint32_t* visited; uint32_t visit_epoch; size_t visited_cap; }; /* ── small helpers ────────────────────────────────────────────────────────── */ static float* vec_normalise_copy(const float* v, int dim){ float* out = (float*)malloc((size_t)dim * sizeof(float)); if (!out) return NULL; double ss = 0.0; for (int i=0;i 0.0){ float inv = (float)(1.0 / sqrt(ss)); for (int i=0;idim; float s0=0,s1=0,s2=0,s3=0; int i=0; for (; i+4<=dim; i+=4){ s0 += a[i]*b[i]; s1 += a[i+1]*b[i+1]; s2 += a[i+2]*b[i+2]; s3 += a[i+3]*b[i+3]; } float dot = (s0+s1)+(s2+s3); for (; icount == nl->cap){ int nc = nl->cap ? nl->cap*2 : 4; int* np = (int*)realloc(nl->ids, (size_t)nc*sizeof(int)); if (!np) return -1; nl->ids = np; nl->cap = nc; } nl->ids[nl->count++] = id; return 0; } /* ── binary heaps over (dist,elem) pairs ──────────────────────────────────── */ typedef struct { float d; int e; } Pair; typedef struct { Pair* a; int n, cap; } Heap; static int heap_reserve(Heap* h, int need){ if (need <= h->cap) return 0; int nc = h->cap ? h->cap*2 : 16; while (nc < need) nc *= 2; Pair* na = (Pair*)realloc(h->a, (size_t)nc*sizeof(Pair)); if (!na) return -1; h->a = na; h->cap = nc; return 0; } /* Order predicate: for a MAX-heap on distance, "higher priority" = larger dist; * ties broken by larger element index (deterministic + stable). is_max selects. */ static inline int pair_before(Pair x, Pair y, int is_max){ if (x.d != y.d) return is_max ? (x.d > y.d) : (x.d < y.d); return is_max ? (x.e > y.e) : (x.e < y.e); } static int heap_push(Heap* h, Pair v, int is_max){ if (heap_reserve(h, h->n+1)) return -1; int i = h->n++; h->a[i] = v; while (i > 0){ int p = (i-1)/2; if (pair_before(h->a[i], h->a[p], is_max)){ Pair t=h->a[i]; h->a[i]=h->a[p]; h->a[p]=t; i=p; } else break; } return 0; } static Pair heap_pop(Heap* h, int is_max){ Pair top = h->a[0]; h->a[0] = h->a[--h->n]; int i = 0; for (;;){ int l=2*i+1, r=2*i+2, best=i; if (ln && pair_before(h->a[l], h->a[best], is_max)) best=l; if (rn && pair_before(h->a[r], h->a[best], is_max)) best=r; if (best==i) break; Pair t=h->a[i]; h->a[i]=h->a[best]; h->a[best]=t; i=best; } return top; } /* ── visited set ──────────────────────────────────────────────────────────── */ static int visited_ensure(VIndex* ix){ if (ix->visited_cap >= ix->cap && ix->visited) return 0; size_t nc = ix->cap ? ix->cap : 16; uint32_t* nv = (uint32_t*)realloc(ix->visited, nc*sizeof(uint32_t)); if (!nv) return -1; if (nc > ix->visited_cap) memset(nv + ix->visited_cap, 0, (nc-ix->visited_cap)*sizeof(uint32_t)); ix->visited = nv; ix->visited_cap = nc; return 0; } static inline void visited_reset(VIndex* ix){ if (++ix->visit_epoch == 0){ /* wrapped: clear all */ memset(ix->visited, 0, ix->visited_cap*sizeof(uint32_t)); ix->visit_epoch = 1; } } static inline int is_visited(VIndex* ix, int e){ return ix->visited[e]==ix->visit_epoch; } static inline void mark_visited(VIndex* ix, int e){ ix->visited[e]=ix->visit_epoch; } /* ── search one layer (Algorithm 2): best-first, ef-bounded ───────────────── */ /* Returns results as an unsorted Heap (max-heap on distance, size<=ef). Caller * owns res->a. `q` is a normalised query. */ static int search_layer(VIndex* ix, const float* q, const int* eps, int neps, int ef, int layer, Heap* res /*out, max-heap*/){ Heap cand = {0,0,0}; /* min-heap: nearest to expand */ res->a=NULL; res->n=0; res->cap=0; visited_reset(ix); for (int i=0;ielems[e].vec); Pair p = { d, e }; if (heap_push(&cand,p,0) || heap_push(res,p,1)){ free(cand.a); return -1; } } while (res->n > ef) heap_pop(res,1); /* trim to ef */ while (cand.n > 0){ Pair c = heap_pop(&cand,0); float worst = res->a[0].d; /* farthest kept result */ if (res->n >= ef && c.d > worst) break; Elem* ce = &ix->elems[c.e]; if (layer <= ce->level){ NeighList* nl = &ce->links[layer]; for (int i=0;icount;i++){ int e = nl->ids[i]; if (is_visited(ix,e)) continue; mark_visited(ix,e); float d = vdist(ix, q, ix->elems[e].vec); if (res->n < ef || d < res->a[0].d){ Pair p = { d, e }; if (heap_push(&cand,p,0) || heap_push(res,p,1)){ free(cand.a); return -1; } if (res->n > ef) heap_pop(res,1); } } } } free(cand.a); return 0; } /* ── neighbour selection heuristic (Algorithm 4) ──────────────────────────── */ /* From candidate pairs W (any order), pick up to M diverse neighbours of q. * Keep c only if it is nearer to q than to every already-chosen neighbour; * backfill from the pruned set (nearest first) to reach M for connectivity. * Writes chosen element indices into out[], returns the count. */ static int select_neighbors(VIndex* ix, const float* q, Pair* W, int nW, int M, int* out){ (void)q; /* q's distances are precomputed in W[].d; kept for call-site clarity */ /* sort W ascending by (dist,elem) — deterministic. */ for (int i=1;i=0 && !pair_before(W[j],key,0)){ W[j+1]=W[j]; j--; } W[j+1]=key; } int nout = 0; Pair* pruned = (Pair*)malloc((size_t)(nW?nW:1)*sizeof(Pair)); int npr = 0; if (!pruned) return -1; for (int i=0;ielems[W[i].e].vec, ix->elems[out[j]].vec); if (d < W[i].d){ good = 0; break; } /* nearer an existing pick → drop */ } if (good) out[nout++] = W[i].e; else pruned[npr++] = W[i]; } for (int i=0;ielems[e].links[layer]; if (nl->count <= Mmax) return; const float* base = ix->elems[e].vec; Pair* W = (Pair*)malloc((size_t)nl->count*sizeof(Pair)); if (!W) return; int nW = nl->count; for (int i=0;ielems[nl->ids[i]].vec), nl->ids[i] }; int* keep = (int*)malloc((size_t)nW*sizeof(int)); if (!keep){ free(W); return; } int nk = select_neighbors(ix, base, W, nW, Mmax, keep); if (nk >= 0){ nl->count = nk; for (int i=0;iids[i]=keep[i]; } free(keep); free(W); } /* ── insert ───────────────────────────────────────────────────────────────── */ static int elems_reserve(VIndex* ix){ if (ix->n < ix->cap) return 0; size_t nc = ix->cap ? ix->cap*2 : 64; Elem* ne = (Elem*)realloc(ix->elems, nc*sizeof(Elem)); if (!ne) return -1; ix->elems = ne; ix->cap = nc; return visited_ensure(ix); } int vindex_insert(VIndex* ix, uint64_t node_id, const float* vec){ if (!ix || !vec) return -1; if (elems_reserve(ix)) return -1; int cur = (int)ix->n; /* deterministic level assignment, seeded per-node. */ uint64_t seed = VINDEX_FIXED_SEED ^ (node_id + 0x2545F4914F6CDD1DULL*(uint64_t)cur); int level = (int)(-log(sm_uniform(&seed)) * ix->mL); if (level < 0) level = 0; Elem* el = &ix->elems[cur]; el->node_id = node_id; el->level = level; el->vec = vec_normalise_copy(vec, ix->dim); el->links = (NeighList*)calloc((size_t)level+1, sizeof(NeighList)); if (!el->vec || !el->links){ free(el->vec); free(el->links); return -1; } ix->n++; if (ix->entry < 0){ /* first element */ ix->entry = cur; ix->max_level = level; return 0; } int ep = ix->entry; int L = ix->max_level; /* greedy descent through layers above `level` to refine the entry point. */ for (int lc = L; lc > level; lc--){ Heap r = {0,0,0}; int eps1[1] = { ep }; if (search_layer(ix, el->vec, eps1, 1, 1, lc, &r)){ return -1; } if (r.n){ ep = r.a[0].e; float bd=r.a[0].d; for (int i=1;i= 0; lc--){ int Mmax = (lc==0) ? ix->M0 : ix->M; Heap W = {0,0,0}; if (search_layer(ix, el->vec, eps, neps, ix->ef_construction, lc, &W)){ rc=-1; break; } int* chosen = (int*)malloc((size_t)(W.n?W.n:1)*sizeof(int)); if (!chosen){ free(W.a); rc=-1; break; } int nc = select_neighbors(ix, el->vec, W.a, W.n, Mmax, chosen); if (nc < 0){ free(chosen); free(W.a); rc=-1; break; } /* link cur <-> chosen (bidirectional), prune neighbours if over-full. */ for (int i=0;ilinks[lc], nb) || nl_push(&ix->elems[nb].links[lc], cur)){ free(chosen); free(W.a); rc=-1; goto done; } prune_links(ix, nb, lc, Mmax); } free(chosen); /* next layer's entry points = this layer's ef results. */ if (lc > 0){ int* neweps = (int*)malloc((size_t)(W.n?W.n:1)*sizeof(int)); if (!neweps){ free(W.a); rc=-1; break; } for (int i=0;i ix->max_level){ ix->max_level = level; ix->entry = cur; } return 0; } /* ── search ───────────────────────────────────────────────────────────────── */ int vindex_search(VIndex* ix, const float* query, int k, int ef_search, uint64_t* node_id_out, float* dist_out){ if (!ix || !query || k <= 0) return -1; if (ix->entry < 0) return 0; if (ef_search <= 0) ef_search = VINDEX_DEFAULT_EF_SEARCH; if (ef_search < k) ef_search = k; float* q = vec_normalise_copy(query, ix->dim); if (!q) return -1; int ep = ix->entry; for (int lc = ix->max_level; lc > 0; lc--){ Heap r = {0,0,0}; int eps[1] = { ep }; if (search_layer(ix, q, eps, 1, 1, lc, &r)){ free(q); return -1; } if (r.n){ int b=r.a[0].e; float bd=r.a[0].d; for (int i=1;i=0;i--) sorted[i] = heap_pop(&res,1); /* farthest first out → fill from end */ free(res.a); int out_n = (k < total) ? k : total; for (int i=0;ielems[sorted[i].e].node_id; if (dist_out) dist_out[i] = sorted[i].d; } free(sorted); return out_n; } size_t vindex_size(const VIndex* ix){ return ix ? ix->n : 0; } VIndex* vindex_create(int dim, int M, int ef_construction){ if (dim <= 0) return NULL; if (M <= 0) M = VINDEX_DEFAULT_M; if (ef_construction <= 0) ef_construction = VINDEX_DEFAULT_EF_CONSTRUCTION; VIndex* ix = (VIndex*)calloc(1, sizeof(VIndex)); if (!ix) return NULL; ix->dim = dim; ix->M = M; ix->M0 = 2*M; ix->ef_construction = ef_construction; ix->mL = 1.0 / log((double)M > 1.0 ? (double)M : 2.0); ix->entry = -1; ix->max_level = 0; ix->visit_epoch = 0; return ix; } void vindex_free(VIndex* ix){ if (!ix) return; for (size_t i=0;in;i++){ Elem* e = &ix->elems[i]; if (e->links) for (int l=0;l<=e->level;l++) free(e->links[l].ids); free(e->links); free(e->vec); } free(ix->elems); free(ix->visited); free(ix); } /* ── read-only decode of the paged store node format (design §2.4) ─────────── */ /* Mirrors engram_store.c constants; the on-disk format is PERMANENT so these are * safe to duplicate for a read-only harvest of emb vectors. */ #define VS_PAGE_SIZE 16384u #define VS_HDR 32u #define VS_SLOT_SIZE 6u #define VS_SLOT_LIVE 1u #define VS_REC_HDR 4u #define VS_REC_OVERFLOW 1u #define VS_PT_NODE 1u #define VS_OVF_NEXT 32u #define VS_OVF_LEN 40u #define VS_OVF_DATA 44u #define VS_NT_ID 1u #define VS_NT_EMB 24u #define VS_NT_EMB_DIM 25u static uint16_t vg_u16(const uint8_t* p){ return (uint16_t)(p[0] | (p[1]<<8)); } static uint32_t vg_u32(const uint8_t* p){ uint32_t v=0; for(int i=0;i<4;i++) v|=(uint32_t)p[i]<<(8*i); return v; } static uint64_t vg_u64(const uint8_t* p){ uint64_t v=0; for(int i=0;i<8;i++) v|=(uint64_t)p[i]<<(8*i); return v; } static int vs_pread(int fd, uint64_t page, uint8_t* buf){ off_t off = (off_t)page * VS_PAGE_SIZE; ssize_t r = pread(fd, buf, VS_PAGE_SIZE, off); return (r == (ssize_t)VS_PAGE_SIZE) ? 0 : -1; } /* Read a (possibly overflowed) record body; caller frees *out. */ static int vs_read_body(int fd, const uint8_t* page, uint16_t off, uint16_t len, uint8_t** out, size_t* outlen){ if (len < VS_REC_HDR) return -1; uint8_t flags = page[off+3]; if (flags & VS_REC_OVERFLOW){ uint64_t head = vg_u64(page + off + VS_REC_HDR); uint64_t total = vg_u64(page + off + VS_REC_HDR + 8); uint8_t* body = (uint8_t*)malloc(total ? total : 1); if (!body) return -1; size_t got=0; uint64_t id=head; uint8_t ov[VS_PAGE_SIZE]; while (id){ if (vs_pread(fd, id, ov)){ free(body); return -1; } uint32_t chunk = vg_u32(ov + VS_OVF_LEN); if (got + chunk > total){ free(body); return -1; } memcpy(body+got, ov+VS_OVF_DATA, chunk); got += chunk; id = vg_u64(ov + VS_OVF_NEXT); } if (got != total){ free(body); return -1; } *out = body; *outlen = total; } else { uint16_t reclen = vg_u16(page + off); if (reclen < VS_REC_HDR) return -1; size_t blen = reclen - VS_REC_HDR; uint8_t* body = (uint8_t*)malloc(blen ? blen : 1); if (!body) return -1; memcpy(body, page + off + VS_REC_HDR, blen); *out = body; *outlen = blen; } return 0; } /* Extract id (strdup) and emb (malloc'd float[dim]) from a TLV node body. */ static void vs_parse_node(const uint8_t* body, size_t len, char** id_out, float** emb_out, int* dim_out){ *id_out=NULL; *emb_out=NULL; *dim_out=0; size_t i=0; while (i + 5 <= len){ uint8_t tag = body[i]; uint32_t flen = vg_u32(body + i + 1); if (i + 5 + (size_t)flen > len) break; const uint8_t* v = body + i + 5; if (tag == VS_NT_ID){ char* s = (char*)malloc(flen+1); if (s){ memcpy(s,v,flen); s[flen]=0; free(*id_out); *id_out=s; } } else if (tag == VS_NT_EMB){ int dim = (int)(flen/4); float* e = (float*)malloc((size_t)(dim?dim:1)*sizeof(float)); if (e){ for (int k=0;kn*2 >= s->cap){ size_t nc = s->cap ? s->cap*2 : 1024; char** nk = (char**)calloc(nc, sizeof(char*)); if (!nk) return -1; for (size_t i=0;icap;i++) if (s->k[i]){ size_t j=vs_fnv(s->k[i])&(nc-1); while(nk[j]) j=(j+1)&(nc-1); nk[j]=s->k[i]; } free(s->k); s->k=nk; s->cap=nc; } size_t j = vs_fnv(key)&(s->cap-1); while (s->k[j]){ if (strcmp(s->k[j],key)==0) return 0; j=(j+1)&(s->cap-1); } char* d = strdup(key); if(!d) return -1; s->k[j]=d; s->n++; return 1; } static void strset_free(StrSet* s){ for(size_t i=0;icap;i++) free(s->k[i]); free(s->k); } int vindex_build_from_store(VIndex* ix, const char* store_path, char*** ids_out, int* n_out){ if (!ix || !store_path) return -1; int fd = open(store_path, O_RDONLY); if (fd < 0) return -1; struct stat st; if (fstat(fd, &st) != 0){ close(fd); return -1; } uint64_t npages = (uint64_t)st.st_size / VS_PAGE_SIZE; char** ids = NULL; size_t ids_n = 0, ids_cap = 0; StrSet seen = {0,0,0}; int inserted = 0; uint8_t page[VS_PAGE_SIZE]; for (uint64_t pg = 2; pg < npages; pg++){ /* pages 0,1 = superblocks */ if (vs_pread(fd, pg, page)) continue; if (page[8] != VS_PT_NODE) continue; int slots = vg_u16(page + 10); for (int sidx=0; sidx VS_PAGE_SIZE) continue; uint8_t* body=NULL; size_t blen=0; if (vs_read_body(fd, page, off, len, &body, &blen)) continue; char* id=NULL; float* emb=NULL; int dim=0; vs_parse_node(body, blen, &id, &emb, &dim); free(body); if (!id || !emb || dim != ix->dim){ free(id); free(emb); continue; } int add = strset_add(&seen, id); if (add <= 0){ free(id); free(emb); continue; } /* dup or err */ if (vindex_insert(ix, (uint64_t)inserted, emb) != 0){ free(id); free(emb); break; } free(emb); if (ids_n == ids_cap){ size_t nc = ids_cap ? ids_cap*2 : 256; char** ni = (char**)realloc(ids, nc*sizeof(char*)); if (!ni){ free(id); break; } ids = ni; ids_cap = nc; } ids[ids_n++] = id; /* transfers ownership */ inserted++; } } close(fd); strset_free(&seen); if (ids_out){ *ids_out = ids; if (n_out) *n_out = (int)ids_n; } else { for (size_t i=0;idim, ix->M, ix->ef_construction, (int32_t)ix->n, ix->entry, ix->max_level }; WR(hdr, sizeof(hdr)); for (size_t i=0; ok && in; i++){ Elem* e = &ix->elems[i]; WR(&e->node_id, sizeof(uint64_t)); int32_t lvl = e->level; WR(&lvl, sizeof(int32_t)); WR(e->vec, (size_t)ix->dim*sizeof(float)); for (int l=0; ok && l<=e->level; l++){ int32_t c = e->links[l].count; WR(&c, sizeof(int32_t)); WR(e->links[l].ids, (size_t)c*sizeof(int)); } } #undef WR fclose(f); return ok ? 0 : -1; } VIndex* vindex_load(const char* path){ FILE* f = fopen(path, "rb"); if (!f) return NULL; char magic[8]; if (fread(magic,1,8,f)!=8 || memcmp(magic,VINDEX_SAVE_MAGIC,8)!=0){ fclose(f); return NULL; } int32_t hdr[6]; if (fread(hdr,sizeof(hdr),1,f)!=1){ fclose(f); return NULL; } VIndex* ix = vindex_create(hdr[0], hdr[1], hdr[2]); if (!ix){ fclose(f); return NULL; } size_t N = (size_t)hdr[3]; int ok = 1; for (size_t i=0; ok && ielems[ix->n]; int32_t lvl; if (fread(&e->node_id,sizeof(uint64_t),1,f)!=1 || fread(&lvl,sizeof(int32_t),1,f)!=1){ ok=0; break; } e->level = lvl; e->vec = (float*)malloc((size_t)ix->dim*sizeof(float)); e->links = (NeighList*)calloc((size_t)lvl+1, sizeof(NeighList)); if (!e->vec || !e->links){ free(e->vec); free(e->links); ok=0; break; } if (fread(e->vec,sizeof(float),(size_t)ix->dim,f)!=(size_t)ix->dim){ ok=0; } for (int l=0; ok && l<=lvl; l++){ int32_t c; if (fread(&c,sizeof(int32_t),1,f)!=1){ ok=0; break; } e->links[l].ids = (int*)malloc((size_t)(c?c:1)*sizeof(int)); e->links[l].cap = c; e->links[l].count = c; if (c && fread(e->links[l].ids,sizeof(int),(size_t)c,f)!=(size_t)c){ ok=0; } } ix->n++; } ix->entry = hdr[4]; ix->max_level = hdr[5]; fclose(f); if (!ok){ vindex_free(ix); return NULL; } return ix; }