From 5336cfe0a6f00cc1b411c3b91218dd8761515f7d Mon Sep 17 00:00:00 2001 From: Will Anderson Date: Thu, 13 Aug 2026 00:26:32 -0500 Subject: [PATCH] =?UTF-8?q?M9=20=C2=A75:=20geometry=20OPERATORS=20as=20C?= =?UTF-8?q?=20functions=20+=20EL=20builtins=20(read-only,=20staged)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bring the relational-neighborhood geometry OPERATORS from the viz proxy (engram-geometry-proxy.py §5) into the C runtime as reusable primitives, and expose each as an EL builtin so any CGI app / el program can use them — not just the engram service internals. engram_geometry.{h,c} (pure, libm-only, read-only over descriptors): - engram_geo_overlap : shared-member Jaccard + centroid/scale proximity score + intersection centroid. - engram_geo_subtract : orthogonal-complement residual (project A onto I - V_B V_Bᵀ), closed-form variance_explained_by_B, residual ellipsoid + centroid-diff; set-diff variant. - engram_geo_combine : pooled descriptor (exact law-of-total-variance mean + covariance), re-eigendecomposed. - engram_geo_distance : centroid L2 + cosine + closed-form Wasserstein-2 (Bures) — mirrors the proxy _wasserstein2. - engram_geo_analogy : orthogonal Procrustes R = UVᵀ (SVD) aligning A's principal frame to B's + apply helper. The C descriptor is full-dim/centered with a low-rank covariance from its top axes; operators mirror the proxy FORMULAS and do the Wasserstein/combine eigen work inside the small joint-axis subspace (exact there). Reuses jacobi_sym. EL builtins (el_runtime.{h,c}, el_seed.c native wrappers): engram_geo_{descriptor,overlap,subtract,combine,distance,analogy}_json — take comma-separated seed-id set(s), build the CENTERED descriptor against the true store-wide mean (ad-hoc path), run the operator, return JSON. Additive: no flag, no effect on activation/retrieval. Surfacing via engram.el + the elc fold is a cutover step (same elc-drift deferral as the P0/P5 builtins); the C table wiring is registered now. Tested: synthetic closed-form unit suite (20/20) — Wasserstein, Jaccard/score, orthogonal residual, set-diff, pooled combine, Procrustes recovery. ASan+UBSan clean; 0 leaks. Read-only; no activation/retrieval behavior change. --- lang/runtime/el_runtime.c | 179 ++++++++++++ lang/runtime/el_runtime.h | 7 + lang/runtime/el_seed.c | 21 ++ lang/runtime/engram_geometry.c | 499 +++++++++++++++++++++++++++++++++ lang/runtime/engram_geometry.h | 102 +++++++ 5 files changed, 808 insertions(+) diff --git a/lang/runtime/el_runtime.c b/lang/runtime/el_runtime.c index 535d3a5..ec1aa20 100644 --- a/lang/runtime/el_runtime.c +++ b/lang/runtime/el_runtime.c @@ -12216,6 +12216,185 @@ el_val_t engram_dreams_json(el_val_t since_ms) { return el_wrap_str(b.buf); } +/* ═══════════════════════════════════════════════════════════════════════════ + * §5 GEOMETRY OPERATORS as EL BUILTINS (Will: "primitives any CGI application + * should be able to use"). READ-ONLY. Each builds the CENTERED descriptor for a + * comma-separated seed-id set (against the true store-wide mean, exactly as the + * ENGRAM_GEO_PRIMING_NOCACHE ad-hoc path does), then runs the pure operator from + * engram_geometry.c and serializes the result to JSON. Additive: no flag, no + * effect on activation/retrieval. Surfacing these through engram.el + the elc fold + * is a CUTOVER step (same elc-drift deferral as the M-INTEROCEPTION P0/P5 builtins); + * the C table wiring (native __ wrappers in el_seed.c) is in place now. + * ═══════════════════════════════════════════════════════════════════════════ */ + +/* Emit a dim-length float vector as a JSON array of %.6g. */ +static void eg_geo_emit_vec(JsonBuf* b, const float* v, int dim) { + char t[48]; jb_putc(b, '['); + for (int i = 0; i < dim; i++) { snprintf(t, sizeof t, "%s%.6g", i ? "," : "", (double)v[i]); jb_puts(b, t); } + jb_putc(b, ']'); +} + +/* Build the centered geometry descriptor for a comma-separated seed-id list. + * Returns a malloc'd descriptor (engram_geo_free) or NULL if geometry is + * unavailable (no paged store / no embeddings / seeds unresolvable). */ +static GeoDescriptor* eg_geo_build_desc(const char* csv) { + if (!csv || !*csv || !g_engram_store) return NULL; + EngramStore* g = engram_get(); + if (!g || g->node_count <= 0) return NULL; + /* split CSV → seed id array (dup, trimmed of spaces). */ + int cap = 8, ns = 0; char** ids = malloc((size_t)cap * sizeof(char*)); + if (!ids) return NULL; + const char* p = csv; + while (*p) { + while (*p == ' ' || *p == ',') p++; + const char* s = p; + while (*p && *p != ',') p++; + const char* e = p; while (e > s && e[-1] == ' ') e--; + if (e > s) { + if (ns == cap) { cap *= 2; char** t = realloc(ids, (size_t)cap * sizeof(char*)); if (!t) break; ids = t; } + ids[ns] = strndup(s, (size_t)(e - s)); ns++; + } + } + if (ns == 0) { free(ids); return NULL; } + /* infer embedding dim from the first resolvable embedded seed. */ + int32_t dim = 0; + for (int i = 0; i < ns && dim == 0; i++) { + int64_t idx = engram_find_node_index(ids[i]); + if (idx >= 0 && idx < g->node_count && g->nodes[idx].emb && g->nodes[idx].emb_dim > 0) + dim = g->nodes[idx].emb_dim; + } + GeoDescriptor* geo = NULL; + const float* gmean = (dim > 0) ? eg_geo_mean_sync(dim) : NULL; + char** vids = malloc((size_t)g->node_count * sizeof(char*)); + if (gmean && vids) { + for (int64_t i = 0; i < g->node_count; i++) vids[i] = g->nodes[i].id; + geo = engram_geometry_descriptor(g_engram_store, _eg_vindex, vids, (int)g->node_count, + (const char* const*)ids, (size_t)ns, NULL, gmean); + } + free(vids); + for (int i = 0; i < ns; i++) free(ids[i]); + free(ids); + return geo; +} + +static el_val_t eg_geo_err(const char* msg) { + JsonBuf b; jb_init(&b); + jb_puts(&b, "{\"error\":"); jb_emit_escaped(&b, msg); jb_putc(&b, '}'); + return el_wrap_str(b.buf); +} + +/* engram_geo_descriptor_json(seed_ids_csv) → the compact centered descriptor. */ +el_val_t engram_geo_descriptor_json(el_val_t seeds) { + GeoDescriptor* g = eg_geo_build_desc(EL_CSTR(seeds)); + if (!g) return eg_geo_err("geometry unavailable (no paged store / embeddings / seeds)"); + JsonBuf b; jb_init(&b); char t[80]; + jb_putc(&b, '{'); + jb_puts(&b, "\"hub_id\":"); jb_emit_escaped(&b, g->hub_id ? g->hub_id : ""); + snprintf(t, sizeof t, ",\"dim\":%d,\"radius\":%.6g,\"total_variance\":%.6g", g->dim, g->radius, g->total_variance); jb_puts(&b, t); + snprintf(t, sizeof t, ",\"k_core\":%d,\"co_registration\":%.6g", g->k_core, g->co_registration); jb_puts(&b, t); + snprintf(t, sizeof t, ",\"n_members\":%d,\"n_embedded\":%d,\"n_axes\":%d", g->n_members, g->n_embedded, g->n_axes); jb_puts(&b, t); + jb_puts(&b, ",\"axis_extents\":["); for (int i = 0; i < g->n_axes; i++) { snprintf(t, sizeof t, "%s%.6g", i ? "," : "", g->axes[i].extent); jb_puts(&b, t); } jb_putc(&b, ']'); + jb_puts(&b, ",\"members\":["); + for (int i = 0; i < g->n_members; i++) { + if (i) jb_putc(&b, ','); + jb_puts(&b, "{\"id\":"); jb_emit_escaped(&b, g->members[i].id ? g->members[i].id : ""); + snprintf(t, sizeof t, ",\"m\":%.4g,\"cent\":%.4g,\"core\":%d}", g->members[i].membership, g->members[i].centrality, g->members[i].core); jb_puts(&b, t); + } + jb_putc(&b, ']'); jb_putc(&b, '}'); + engram_geo_free(g); + return el_wrap_str(b.buf); +} + +/* engram_geo_overlap_json(a_csv, b_csv) → shared set + Jaccard + score. */ +el_val_t engram_geo_overlap_json(el_val_t a_seeds, el_val_t b_seeds) { + GeoDescriptor* A = eg_geo_build_desc(EL_CSTR(a_seeds)); + GeoDescriptor* B = eg_geo_build_desc(EL_CSTR(b_seeds)); + if (!A || !B) { if (A) engram_geo_free(A); if (B) engram_geo_free(B); return eg_geo_err("geometry unavailable"); } + GeoOverlap o; char t[80]; JsonBuf b; jb_init(&b); + if (engram_geo_overlap(A, B, &o) != 0) { engram_geo_free(A); engram_geo_free(B); return eg_geo_err("dim mismatch"); } + jb_putc(&b, '{'); + snprintf(t, sizeof t, "\"n_shared\":%d,\"n_union\":%d,\"jaccard\":%.6g", o.n_shared, o.n_union, o.jaccard); jb_puts(&b, t); + snprintf(t, sizeof t, ",\"centroid_distance\":%.6g,\"overlap_score\":%.6g", o.centroid_distance, o.overlap_score); jb_puts(&b, t); + jb_puts(&b, ",\"shared_ids\":["); for (int i = 0; i < o.n_shared; i++) { if (i) jb_putc(&b, ','); jb_emit_escaped(&b, o.shared_ids[i]); } jb_putc(&b, ']'); + jb_putc(&b, '}'); + engram_geo_overlap_free(&o); engram_geo_free(A); engram_geo_free(B); + return el_wrap_str(b.buf); +} + +/* engram_geo_subtract_json(a_csv, b_csv, mode) — mode "setdiff" → set-difference, + * anything else → orthogonal-complement residual. */ +el_val_t engram_geo_subtract_json(el_val_t a_seeds, el_val_t b_seeds, el_val_t mode) { + GeoDescriptor* A = eg_geo_build_desc(EL_CSTR(a_seeds)); + GeoDescriptor* B = eg_geo_build_desc(EL_CSTR(b_seeds)); + if (!A || !B) { if (A) engram_geo_free(A); if (B) engram_geo_free(B); return eg_geo_err("geometry unavailable"); } + const char* m = EL_CSTR(mode); char t[80]; JsonBuf b; jb_init(&b); + if (m && strcmp(m, "setdiff") == 0) { + GeoSetDiff s; + if (engram_geo_setdiff(A, B, &s) != 0) { engram_geo_free(A); engram_geo_free(B); return eg_geo_err("dim mismatch"); } + jb_putc(&b, '{'); jb_puts(&b, "\"mode\":\"setdiff\""); + snprintf(t, sizeof t, ",\"n_only\":%d,\"removed\":%d,\"centroid_diff_mag\":%.6g", s.n_only, s.removed, s.centroid_diff_mag); jb_puts(&b, t); + jb_puts(&b, ",\"only_ids\":["); for (int i = 0; i < s.n_only; i++) { if (i) jb_putc(&b, ','); jb_emit_escaped(&b, s.only_ids[i]); } jb_putc(&b, ']'); + jb_putc(&b, '}'); + engram_geo_setdiff_free(&s); + } else { + GeoResidual r; + if (engram_geo_subtract(A, B, 0, &r) != 0) { engram_geo_free(A); engram_geo_free(B); return eg_geo_err("dim mismatch"); } + jb_putc(&b, '{'); jb_puts(&b, "\"mode\":\"residual\""); + snprintf(t, sizeof t, ",\"variance_explained_by_B\":%.6g,\"residual_scale\":%.6g", r.variance_explained_by_B, r.residual_scale); jb_puts(&b, t); + snprintf(t, sizeof t, ",\"removed_dims\":%d,\"residual_n_axes\":%d,\"centroid_diff_mag\":%.6g", r.removed_dims, r.n_axes, r.centroid_diff_mag); jb_puts(&b, t); + jb_putc(&b, '}'); + engram_geo_residual_free(&r); + } + engram_geo_free(A); engram_geo_free(B); + return el_wrap_str(b.buf); +} + +/* engram_geo_combine_json(a_csv, b_csv) → merged descriptor summary. */ +el_val_t engram_geo_combine_json(el_val_t a_seeds, el_val_t b_seeds) { + GeoDescriptor* A = eg_geo_build_desc(EL_CSTR(a_seeds)); + GeoDescriptor* B = eg_geo_build_desc(EL_CSTR(b_seeds)); + if (!A || !B) { if (A) engram_geo_free(A); if (B) engram_geo_free(B); return eg_geo_err("geometry unavailable"); } + GeoDescriptor* C = engram_geo_combine(A, B, 8); + engram_geo_free(A); engram_geo_free(B); + if (!C) return eg_geo_err("combine failed (dim mismatch / OOM)"); + JsonBuf b; jb_init(&b); char t[80]; + jb_putc(&b, '{'); + jb_puts(&b, "\"hub_id\":"); jb_emit_escaped(&b, C->hub_id ? C->hub_id : ""); + snprintf(t, sizeof t, ",\"dim\":%d,\"radius\":%.6g,\"total_variance\":%.6g", C->dim, C->radius, C->total_variance); jb_puts(&b, t); + snprintf(t, sizeof t, ",\"n_members\":%d,\"n_embedded\":%d,\"n_axes\":%d", C->n_members, C->n_embedded, C->n_axes); jb_puts(&b, t); + jb_puts(&b, ",\"axis_extents\":["); for (int i = 0; i < C->n_axes; i++) { snprintf(t, sizeof t, "%s%.6g", i ? "," : "", C->axes[i].extent); jb_puts(&b, t); } jb_putc(&b, ']'); + jb_putc(&b, '}'); + engram_geo_free(C); + return el_wrap_str(b.buf); +} + +/* engram_geo_distance_json(a_csv, b_csv) → centroid + Wasserstein-2. */ +el_val_t engram_geo_distance_json(el_val_t a_seeds, el_val_t b_seeds) { + GeoDescriptor* A = eg_geo_build_desc(EL_CSTR(a_seeds)); + GeoDescriptor* B = eg_geo_build_desc(EL_CSTR(b_seeds)); + if (!A || !B) { if (A) engram_geo_free(A); if (B) engram_geo_free(B); return eg_geo_err("geometry unavailable"); } + GeoDistance d; char t[96]; JsonBuf b; jb_init(&b); + if (engram_geo_distance(A, B, &d) != 0) { engram_geo_free(A); engram_geo_free(B); return eg_geo_err("dim mismatch"); } + snprintf(t, sizeof t, "{\"centroid_distance\":%.6g,\"centroid_cosine\":%.6g,\"wasserstein2\":%.6g}", d.centroid_distance, d.centroid_cosine, d.wasserstein2); + jb_puts(&b, t); + engram_geo_free(A); engram_geo_free(B); + return el_wrap_str(b.buf); +} + +/* engram_geo_analogy_json(a_csv, b_csv) → orthogonal Procrustes residual + rank. */ +el_val_t engram_geo_analogy_json(el_val_t a_seeds, el_val_t b_seeds) { + GeoDescriptor* A = eg_geo_build_desc(EL_CSTR(a_seeds)); + GeoDescriptor* B = eg_geo_build_desc(EL_CSTR(b_seeds)); + if (!A || !B) { if (A) engram_geo_free(A); if (B) engram_geo_free(B); return eg_geo_err("geometry unavailable"); } + GeoAnalogy an; char t[96]; JsonBuf b; jb_init(&b); + if (engram_geo_analogy(A, B, &an) != 0) { engram_geo_free(A); engram_geo_free(B); return eg_geo_err("dim mismatch"); } + snprintf(t, sizeof t, "{\"subspace_rank\":%d,\"residual\":%.6g}", an.r, an.residual); + jb_puts(&b, t); + engram_geo_analogy_free(&an); + engram_geo_free(A); engram_geo_free(B); + return el_wrap_str(b.buf); +} + el_val_t engram_neighbors_json(el_val_t node_id, el_val_t max_depth, el_val_t direction) { /* Re-implement here directly so we serialize without going through * the ElList path. Walks BFS to max_depth, emits {node, edge, hops} diff --git a/lang/runtime/el_runtime.h b/lang/runtime/el_runtime.h index 73ad589..0dd4ef8 100644 --- a/lang/runtime/el_runtime.h +++ b/lang/runtime/el_runtime.h @@ -623,6 +623,13 @@ el_val_t engram_scan_nodes_json(el_val_t limit, el_val_t offset); el_val_t engram_scan_nodes_by_type_json(el_val_t node_type, el_val_t limit, el_val_t offset); el_val_t engram_scan_nodes_emb_json(el_val_t limit, el_val_t offset); el_val_t engram_dreams_json(el_val_t since_ms); +/* §5 geometry operators as EL builtins (read-only; seed-id CSV args). */ +el_val_t engram_geo_descriptor_json(el_val_t seeds); +el_val_t engram_geo_overlap_json(el_val_t a_seeds, el_val_t b_seeds); +el_val_t engram_geo_subtract_json(el_val_t a_seeds, el_val_t b_seeds, el_val_t mode); +el_val_t engram_geo_combine_json(el_val_t a_seeds, el_val_t b_seeds); +el_val_t engram_geo_distance_json(el_val_t a_seeds, el_val_t b_seeds); +el_val_t engram_geo_analogy_json(el_val_t a_seeds, el_val_t b_seeds); el_val_t engram_consolidate_permanence(el_val_t node_id); el_val_t engram_age_field(el_val_t delta_ms); el_val_t engram_age_field_catchup(void); diff --git a/lang/runtime/el_seed.c b/lang/runtime/el_seed.c index 9dfa88b..0a96e00 100644 --- a/lang/runtime/el_seed.c +++ b/lang/runtime/el_seed.c @@ -1094,6 +1094,27 @@ el_val_t __engram_dreams_json(el_val_t since_ms) { return engram_dreams_json(since_ms); } +/* §5 geometry operators — native wrappers (surfacing via engram.el + elc fold is + * the cutover step; the C table wiring is registered here now, per P0/P5). */ +el_val_t __engram_geo_descriptor_json(el_val_t seeds) { + return engram_geo_descriptor_json(seeds); +} +el_val_t __engram_geo_overlap_json(el_val_t a_seeds, el_val_t b_seeds) { + return engram_geo_overlap_json(a_seeds, b_seeds); +} +el_val_t __engram_geo_subtract_json(el_val_t a_seeds, el_val_t b_seeds, el_val_t mode) { + return engram_geo_subtract_json(a_seeds, b_seeds, mode); +} +el_val_t __engram_geo_combine_json(el_val_t a_seeds, el_val_t b_seeds) { + return engram_geo_combine_json(a_seeds, b_seeds); +} +el_val_t __engram_geo_distance_json(el_val_t a_seeds, el_val_t b_seeds) { + return engram_geo_distance_json(a_seeds, b_seeds); +} +el_val_t __engram_geo_analogy_json(el_val_t a_seeds, el_val_t b_seeds) { + return engram_geo_analogy_json(a_seeds, b_seeds); +} + el_val_t __engram_consolidate_permanence(el_val_t node_id) { return engram_consolidate_permanence(node_id); } diff --git a/lang/runtime/engram_geometry.c b/lang/runtime/engram_geometry.c index 5efbb3a..c559939 100644 --- a/lang/runtime/engram_geometry.c +++ b/lang/runtime/engram_geometry.c @@ -585,6 +585,505 @@ void engram_geo_displacement(const GeoDescriptor* a, const GeoDescriptor* b, out->periph_disp = periph_n ? periph_sum/periph_n : 0.0; } +/* ═══════════════════════════════════════════════════════════════════════════ + * §5 GEOMETRY OPERATORS — relational algebra over descriptors. Pure, read-only. + * See engram_geometry.h for the frame contract + the low-rank representation note. + * All eigen-work reuses the file-static jacobi_sym above. + * ═══════════════════════════════════════════════════════════════════════════ */ + +/* dot of two dim-length float vectors. */ +static double geo_vdot(const float* a, const float* b, int dim){ + double s=0; for(int d=0;daxis[d]*w[d]; return s; +} +/* Σ_g · w → out (dim). Σ_g = Σ_k extent_k² a_k a_kᵀ (low-rank from the axes). */ +static void geo_cov_apply(const GeoDescriptor* g, const float* w, float* out, int dim){ + for(int d=0;dn_axes;k++){ + double e=g->axes[k].extent; + double c=geo_axis_coef(&g->axes[k],w,dim)*e*e; + const float* a=g->axes[k].axis; + for(int d=0;dn_axes;k++) s+=g->axes[k].extent*g->axes[k].extent; return s; +} +/* Orthonormal basis (modified Gram–Schmidt) of span(cand[0..nc-1]); each cand is + * dim floats. Writes up to nc rows into Q (row-major, caller allocs nc*dim floats). + * Returns the rank r ≤ nc (near-dependent vectors are dropped). */ +static int geo_orthobasis(float* const* cand, int nc, int dim, float* Q){ + int r=0; + double* v=malloc((size_t)dim*sizeof(double)); + if(!v) return 0; + for(int i=0;i1e-6){ for(int d=0;d0?sqrt(w[k]):0.0; s+=V[i*n+k]*sq*V[j*n+k]; } + R[i*n+j]=s; + } + free(A); free(w); free(V); +} + +/* ── overlap ─────────────────────────────────────────────────────────────── */ +int engram_geo_overlap(const GeoDescriptor* a, const GeoDescriptor* b, GeoOverlap* out){ + if(!a||!b||!out||a->dim!=b->dim) return -1; + memset(out,0,sizeof*out); + int dim=a->dim; out->dim=dim; + int cap = a->n_membersn_members ? a->n_members : b->n_members; + out->shared_ids = cap? calloc((size_t)cap,sizeof(char*)) : NULL; + int ns=0; + for(int i=0;in_members;i++){ + const char* id=a->members[i].id; if(!id) continue; + for(int j=0;jn_members;j++){ + if(b->members[j].id && strcmp(b->members[j].id,id)==0){ + if(out->shared_ids) out->shared_ids[ns]=strdup(id); + ns++; break; + } + } + } + out->n_shared=ns; + out->n_union = a->n_members + b->n_members - ns; + out->jaccard = out->n_union? (double)ns/(double)out->n_union : 0.0; + double d=0; + if(a->centroid && b->centroid){ + for(int k=0;kcentroid[k]-b->centroid[k]; d+=df*df; } + d=sqrt(d); + out->intersection_centroid=malloc((size_t)dim*sizeof(float)); + if(out->intersection_centroid) + for(int k=0;kintersection_centroid[k]=0.5f*(a->centroid[k]+b->centroid[k]); + } + out->centroid_distance=d; + double denom=a->radius+b->radius+1e-9; + double prox=1.0 - d/denom; if(prox<0) prox=0; + out->overlap_score = out->jaccard*0.5 + prox*0.5; + return 0; +} +void engram_geo_overlap_free(GeoOverlap* o){ + if(!o) return; + for(int i=0;in_shared;i++) free(o->shared_ids[i]); + free(o->shared_ids); free(o->intersection_centroid); + memset(o,0,sizeof*o); +} + +/* ── subtract: orthogonal-complement residual ────────────────────────────── */ +int engram_geo_subtract(const GeoDescriptor* a, const GeoDescriptor* b, + int b_dims, GeoResidual* out){ + if(!a||!b||!out||a->dim!=b->dim) return -1; + memset(out,0,sizeof*out); + int dim=a->dim; out->dim=dim; + int mB = (b_dims>0) ? b_dims : (b->n_axes<3 ? b->n_axes : 3); + if(mB>b->n_axes) mB=b->n_axes; if(mB<0) mB=0; + out->removed_dims=mB; + + double cA_norm2 = a->centroid ? geo_vdot(a->centroid,a->centroid,dim) : 0.0; + double Qc2=0; + out->residual_centroid = a->centroid ? malloc((size_t)dim*sizeof(float)) : NULL; + if(a->centroid && out->residual_centroid){ + memcpy(out->residual_centroid,a->centroid,(size_t)dim*sizeof(float)); + for(int k=0;kaxes[k],a->centroid,dim); + Qc2 += coef*coef; + const float* ax=b->axes[k].axis; + for(int d=0;dresidual_centroid[d]-=(float)(coef*ax[d]); + } + } + /* variance of A explained by B: Tr(QΣ_A)=Σ_k Σ_j eA_j²(a_k·u_j)². */ + double trA=geo_cov_trace(a), trQA=0; + for(int k=0;kn_axes;j++){ + double c=geo_axis_coef(&b->axes[k],a->axes[j].axis,dim); + double e=a->axes[j].extent; + trQA += e*e*c*c; + } + } + double Etot=cA_norm2+trA, Eres=(cA_norm2-Qc2)+(trA-trQA); + double ex = Etot>1e-12 ? 1.0-Eres/Etot : 0.0; + if(ex<0)ex=0; if(ex>1)ex=1; + out->variance_explained_by_B=ex; + double resvar=trA-trQA; if(resvar<0) resvar=0; + out->residual_scale=sqrt(resvar); + /* residual axes = P⊥ u_j (projected out of B's subspace). */ + if(a->n_axes>0){ + out->axes=calloc((size_t)a->n_axes,sizeof(GeoAxis)); + int na=0; + for(int j=0;jn_axes && out->axes;j++){ + float* r=malloc((size_t)dim*sizeof(float)); if(!r) break; + memcpy(r,a->axes[j].axis,(size_t)dim*sizeof(float)); + for(int k=0;kaxes[k],a->axes[j].axis,dim); + const float* ax=b->axes[k].axis; + for(int d=0;d1e-6){ for(int d=0;daxes[na].axis=r; out->axes[na].extent=a->axes[j].extent*nrm; na++; } + else free(r); + } + out->n_axes=na; + if(na==0){ free(out->axes); out->axes=NULL; } + } + out->centroid_diff = malloc((size_t)dim*sizeof(float)); + double mag=0; + if(a->centroid && b->centroid && out->centroid_diff) + for(int d=0;dcentroid[d]-b->centroid[d]; out->centroid_diff[d]=df; mag+=(double)df*df; } + out->centroid_diff_mag=sqrt(mag); + return 0; +} +void engram_geo_residual_free(GeoResidual* r){ + if(!r) return; + free(r->residual_centroid); free(r->centroid_diff); + for(int i=0;in_axes;i++) free(r->axes[i].axis); + free(r->axes); + memset(r,0,sizeof*r); +} + +/* ── set-difference variant ──────────────────────────────────────────────── */ +int engram_geo_setdiff(const GeoDescriptor* a, const GeoDescriptor* b, GeoSetDiff* out){ + if(!a||!b||!out||a->dim!=b->dim) return -1; + memset(out,0,sizeof*out); + int dim=a->dim; out->dim=dim; + out->only_ids = a->n_members? calloc((size_t)a->n_members,sizeof(char*)) : NULL; + int no=0, rem=0; + for(int i=0;in_members;i++){ + const char* id=a->members[i].id; if(!id) continue; + int in=0; + for(int j=0;jn_members;j++) if(b->members[j].id && strcmp(b->members[j].id,id)==0){ in=1; break; } + if(in) rem++; + else if(out->only_ids){ out->only_ids[no++]=strdup(id); } + else no++; + } + out->n_only=no; out->removed=rem; + out->centroid_diff=malloc((size_t)dim*sizeof(float)); + double mag=0; + if(a->centroid && b->centroid && out->centroid_diff) + for(int d=0;dcentroid[d]-b->centroid[d]; out->centroid_diff[d]=df; mag+=(double)df*df; } + out->centroid_diff_mag=sqrt(mag); + return 0; +} +void engram_geo_setdiff_free(GeoSetDiff* s){ + if(!s) return; + for(int i=0;in_only;i++) free(s->only_ids[i]); + free(s->only_ids); free(s->centroid_diff); + memset(s,0,sizeof*s); +} + +/* ── combine: pooled Gaussian (exact law-of-total-variance) ───────────────── */ +GeoDescriptor* engram_geo_combine(const GeoDescriptor* a, const GeoDescriptor* b, int top_axes){ + if(!a||!b||a->dim!=b->dim) return NULL; + int dim=a->dim; + if(top_axes<=0) top_axes=8; + double nA=a->n_embedded>0?a->n_embedded:a->n_members; + double nB=b->n_embedded>0?b->n_embedded:b->n_members; + if(nA<1) nA=1; if(nB<1) nB=1; + double nt=nA+nB, wA=nA/nt, wB=nB/nt, cross=nA*nB/(nt*nt); + + GeoDescriptor* g=calloc(1,sizeof(GeoDescriptor)); + if(!g) return NULL; + g->dim=dim; + /* pooled centroid (both must be embedded to have a meaningful centroid) */ + float* d=NULL; double dnorm2=0; + if(a->centroid && b->centroid){ + g->centroid=malloc((size_t)dim*sizeof(float)); + d=malloc((size_t)dim*sizeof(float)); + if(!g->centroid||!d){ free(d); engram_geo_free(g); return NULL; } + for(int i=0;icentroid[i]=(float)(wA*a->centroid[i]+wB*b->centroid[i]); + d[i]=a->centroid[i]-b->centroid[i]; dnorm2+=(double)d[i]*d[i]; + } + } + if(a->global_mean){ + g->global_mean=malloc((size_t)dim*sizeof(float)); + if(g->global_mean) memcpy(g->global_mean,a->global_mean,(size_t)dim*sizeof(float)); + } + /* pooled total variance = wA·trA + wB·trB + cross·‖d‖² */ + double trA=geo_cov_trace(a), trB=geo_cov_trace(b); + g->total_variance = wA*trA + wB*trB + cross*dnorm2; + g->radius = sqrt(g->total_variance>0?g->total_variance:0); + + /* eigendecompose the pooled covariance inside the joint subspace. */ + int nc=a->n_axes+b->n_axes+(d?1:0); + if(nc>0){ + float** cand=malloc((size_t)nc*sizeof(float*)); int ci=0; + for(int k=0;kn_axes;k++) cand[ci++]=a->axes[k].axis; + for(int k=0;kn_axes;k++) cand[ci++]=b->axes[k].axis; + if(d) cand[ci++]=d; + float* Q=malloc((size_t)nc*dim*sizeof(float)); + int r=(cand&&Q)?geo_orthobasis(cand,nc,dim,Q):0; + if(r>0){ + /* M[i][j] = q_iᵀ Σ_pooled q_j */ + double* M=calloc((size_t)r*r,sizeof(double)); + float* sqa=malloc((size_t)dim*sizeof(float)); + float* sqb=malloc((size_t)dim*sizeof(float)); + if(M&&sqa&&sqb){ + for(int j=0;jw[ord[i]]){int t=ord[i];ord[i]=ord[j];ord[j]=t;} + int keep=top_axes; if(keep>r) keep=r; + g->axes=calloc((size_t)keep,sizeof(GeoAxis)); + int na=0; + for(int t=0;taxes;t++){ + int c=ord[t]; double lam=w[c]; if(lam<0) lam=0; + if(lam<1e-12) continue; + float* ax=calloc((size_t)dim,sizeof(float)); if(!ax) break; + for(int dd=0;dd1e-12) for(int dd=0;ddaxes[na].axis=ax; g->axes[na].extent=sqrt(lam); na++; + } + g->n_axes=na; + free(ord); + } + free(w); free(V); + } + free(M); free(sqa); free(sqb); + } + free(Q); free(cand); + } + free(d); + + /* member id-union (membership = max of the two copies). */ + int cap=a->n_members+b->n_members; + g->members = cap? calloc((size_t)cap,sizeof(GeoMember)) : NULL; + int M=0; + for(int i=0;in_members && g->members;i++){ + g->members[M].id=strdup(a->members[i].id?a->members[i].id:""); + g->members[M].membership=a->members[i].membership; + g->members[M].centrality=a->members[i].centrality; + g->members[M].salience=a->members[i].salience; + g->members[M].core=a->members[i].core; + g->members[M].dist_centroid=a->members[i].dist_centroid; + g->members[M].embedded=a->members[i].embedded; + M++; + } + for(int j=0;jn_members && g->members;j++){ + const char* id=b->members[j].id; int found=-1; + for(int i=0;imembers[i].id && id && strcmp(g->members[i].id,id)==0){ found=i; break; } + if(found>=0){ + if(b->members[j].membership>g->members[found].membership) + g->members[found].membership=b->members[j].membership; + if(b->members[j].centrality>g->members[found].centrality) + g->members[found].centrality=b->members[j].centrality; + } else { + g->members[M].id=strdup(id?id:""); + g->members[M].membership=b->members[j].membership; + g->members[M].centrality=b->members[j].centrality; + g->members[M].salience=b->members[j].salience; + g->members[M].core=b->members[j].core; + g->members[M].dist_centroid=b->members[j].dist_centroid; + g->members[M].embedded=b->members[j].embedded; + M++; + } + } + g->n_members=M; + g->n_embedded=a->n_embedded+b->n_embedded; + /* hub = highest-centrality union member (fallback A's hub). */ + int hub=-1; double hv=-1; + for(int i=0;imembers[i].centrality>hv){ hv=g->members[i].centrality; hub=i; } + g->hub_id = strdup(hub>=0 ? g->members[hub].id : (a->hub_id?a->hub_id:"")); + g->k_core = a->k_core>b->k_core ? a->k_core : b->k_core; + g->co_registration = 0.5*(a->co_registration+b->co_registration); + g->n_edges=0; g->edges=NULL; + return g; +} + +/* ── distance: centroid + Wasserstein-2 (Bures) ──────────────────────────── */ +int engram_geo_distance(const GeoDescriptor* a, const GeoDescriptor* b, GeoDistance* out){ + if(!a||!b||!out||a->dim!=b->dim) return -1; + memset(out,0,sizeof*out); + int dim=a->dim; out->dim=dim; + double d2=0, cdot=0, na=0, nb=0; + if(a->centroid && b->centroid){ + for(int k=0;kcentroid[k], y=b->centroid[k]; + double df=x-y; d2+=df*df; cdot+=x*y; na+=x*x; nb+=y*y; } + } + out->centroid_distance=sqrt(d2); + out->centroid_cosine = (na>1e-12 && nb>1e-12) ? cdot/(sqrt(na)*sqrt(nb)) : 0.0; + + double trace_term=0; + int nc=a->n_axes+b->n_axes; + if(nc>0){ + float** cand=malloc((size_t)nc*sizeof(float*)); int ci=0; + for(int k=0;kn_axes;k++) cand[ci++]=a->axes[k].axis; + for(int k=0;kn_axes;k++) cand[ci++]=b->axes[k].axis; + float* Q=malloc((size_t)nc*dim*sizeof(float)); + int r=(cand&&Q)?geo_orthobasis(cand,nc,dim,Q):0; + if(r>0){ + double* C1=malloc((size_t)r*r*sizeof(double)); + double* C2=malloc((size_t)r*r*sizeof(double)); + double* s2=malloc((size_t)r*r*sizeof(double)); + double* tmp=malloc((size_t)r*r*sizeof(double)); + double* mid=malloc((size_t)r*r*sizeof(double)); + double* inner=malloc((size_t)r*r*sizeof(double)); + if(C1&&C2&&s2&&tmp&&mid&&inner){ + geo_cov_in_basis(a,Q,r,dim,C1); + geo_cov_in_basis(b,Q,r,dim,C2); + geo_symsqrt(C2,s2,r); /* s2 = C2^{1/2} */ + geo_matmul(s2,C1,tmp,r); geo_matmul(tmp,s2,mid,r); /* s2 C1 s2 */ + geo_symsqrt(mid,inner,r); /* inner = (s2 C1 s2)^{1/2}*/ + double trc=0; + for(int i=0;iwasserstein2=sqrt(w2); + return 0; +} + +/* ── analogy: orthogonal Procrustes (SVD via jacobi on MᵀM) ───────────────── */ +int engram_geo_analogy(const GeoDescriptor* a, const GeoDescriptor* b, GeoAnalogy* out){ + if(!a||!b||!out||a->dim!=b->dim) return -1; + memset(out,0,sizeof*out); + int dim=a->dim; out->dim=dim; + int k = a->n_axesn_axes ? a->n_axes : b->n_axes; /* paired axes */ + if(k<=0){ out->r=0; out->residual=0; return 0; } + int nc=a->n_axes+b->n_axes; + float** cand=malloc((size_t)nc*sizeof(float*)); int ci=0; + for(int t=0;tn_axes;t++) cand[ci++]=a->axes[t].axis; + for(int t=0;tn_axes;t++) cand[ci++]=b->axes[t].axis; + float* Q=malloc((size_t)nc*dim*sizeof(float)); + int r=(cand&&Q)?geo_orthobasis(cand,nc,dim,Q):0; + if(r<=0){ free(cand); free(Q); out->r=0; return 0; } + /* extent-scaled frame coords in Q: Ahat,Bhat are r×k. */ + double* Ah=calloc((size_t)r*k,sizeof(double)); + double* Bh=calloc((size_t)r*k,sizeof(double)); + for(int c=0;caxes[c].extent, eb=b->axes[c].extent; + for(int i=0;iaxes[c],&Q[(size_t)i*dim],dim); + Bh[i*k+c]=eb*geo_axis_coef(&b->axes[c],&Q[(size_t)i*dim],dim); + } + } + /* Mhat = Ah·Bhᵀ (r×r) */ + double* Mh=calloc((size_t)r*r,sizeof(double)); + for(int i=0;i0?sqrt(w[c]):0.0; + if(sig>1e-9){ + for(int i=0;iR=calloc((size_t)r*r,sizeof(double)); + for(int i=0;iR[i*r+j]=s; } + /* residual = ‖Ah − R̂·Bh‖_F */ + double resid=0; + for(int c=0;cR[i*r+j]*Bh[j*k+c]; + double df=Ah[i*k+c]-rb; resid+=df*df; + } + out->residual=sqrt(resid); + out->r=r; + out->basis=malloc((size_t)r*dim*sizeof(float)); + if(out->basis) memcpy(out->basis,Q,(size_t)r*dim*sizeof(float)); + free(Ah); free(Bh); free(Mh); free(S); free(w); free(V); free(U); + free(Q); free(cand); + return 0; +} +void engram_geo_analogy_apply(const GeoAnalogy* an, const float* v, float* out_vec){ + if(!an||!v||!out_vec) return; + int dim=an->dim, r=an->r; + for(int d=0;dbasis||!an->R) return; + double* c=malloc((size_t)r*sizeof(double)); + double* cp=malloc((size_t)r*sizeof(double)); + if(!c||!cp){ free(c); free(cp); return; } + for(int i=0;ibasis[(size_t)i*dim]; + for(int d=0;dR[i*r+j]*c[j]; cp[i]=s; } + for(int i=0;ibasis[(size_t)i*dim]; + for(int d=0;dbasis); free(an->R); + memset(an,0,sizeof*an); +} + /* ═══════════════════════════════════════════════════════════════════════════ * M10 — REIFICATION: persist / load / lookup first-class neighborhood records. * ═══════════════════════════════════════════════════════════════════════════ */ diff --git a/lang/runtime/engram_geometry.h b/lang/runtime/engram_geometry.h index 94fc2e8..00e7110 100644 --- a/lang/runtime/engram_geometry.h +++ b/lang/runtime/engram_geometry.h @@ -174,6 +174,108 @@ typedef struct { void engram_geo_displacement(const GeoDescriptor* a, const GeoDescriptor* b, double core_frac, GeoDisplacement* out); +/* ═══════════════════════════════════════════════════════════════════════════ + * §5 GEOMETRY OPERATORS — a relational ALGEBRA over neighborhood descriptors. + * These are the reusable primitives Will specified: "primitives any CGI + * application should be able to use." READ-ONLY and PURE (stdlib + libm only) — + * they consume GeoDescriptor(s) and never touch the store, index, or activation. + * + * FRAME CONTRACT: both inputs MUST have been built in the SAME frame — identical + * emb `dim` and identical `global_mean` (centered against the one true store-wide + * mean). The reify path builds every neighborhood that way, so descriptors are + * directly comparable. An operator returns <0 / NULL if the dims disagree. + * + * REPRESENTATION: the C descriptor lives in the FULL emb dim with a LOW-RANK + * covariance Σ = Σ_k extent_k² · a_k a_kᵀ over its retained principal axes + * (top_axes; the discarded tail variance is not modeled). Every operator mirrors + * the viz-proxy (engram-geometry-proxy.py §5) FORMULA exactly, but evaluates it on + * this representation — so semantics match the proxy while absolute numbers differ + * (proxy works in a 24-dim global-PCA reduced dense frame; C in full-dim low-rank). + * The Wasserstein / combine eigen-work is done inside the small JOINT axis subspace + * (dimension ≤ nA+nB+1), which is EXACT for the low-rank covariances there. + * Each result struct is released by its engram_geo_*_free. + * ═══════════════════════════════════════════════════════════════════════════ */ + +/* overlap(A,B): shared-member set + Jaccard + centroid/scale proximity score. */ +typedef struct { + char** shared_ids; /* ids present in BOTH neighborhoods (owned) */ + int n_shared; + int n_union; /* |A ∪ B| by id */ + double jaccard; /* |A∩B| / |A∪B| */ + double centroid_distance; /* L2 between the (centered) centroids */ + double overlap_score; /* jacc*0.5 + max(0,1−d/(rA+rB))*0.5 (proxy form)*/ + float* intersection_centroid; /* midpoint of the two centroids (dim, owned) */ + int dim; +} GeoOverlap; +int engram_geo_overlap(const GeoDescriptor* a, const GeoDescriptor* b, GeoOverlap* out); +void engram_geo_overlap_free(GeoOverlap* o); + +/* subtract(A,B) — ORTHOGONAL-COMPLEMENT residual: project A onto I − V_B V_Bᵀ + * (V_B = B's top `b_dims` principal axes) — "A with B's framing removed". Returns + * A's residual centroid + residual ellipsoid, the fraction of A's energy that lives + * inside B's subspace, and the centroid-difference vector. b_dims<=0 → min(3,nB). */ +typedef struct { + int dim; + float* residual_centroid; /* P⊥ c_A (owned) */ + float* centroid_diff; /* c_A − c_B (owned) */ + double centroid_diff_mag; + double variance_explained_by_B; /* (‖Qc_A‖²+Tr(QΣ_A)) / (‖c_A‖²+Tr Σ_A) ∈[0,1]*/ + int removed_dims; /* # of B axes used as V_B */ + double residual_scale; /* sqrt(Tr(P⊥ Σ_A P⊥)) */ + int n_axes; /* residual principal axes (owned) */ + GeoAxis* axes; +} GeoResidual; +int engram_geo_subtract(const GeoDescriptor* a, const GeoDescriptor* b, + int b_dims, GeoResidual* out); +void engram_geo_residual_free(GeoResidual* r); + +/* set-diff variant of subtract: members in A but not in B + the centroid arrow. */ +typedef struct { + char** only_ids; /* member ids in A and not in B (owned) */ + int n_only; + int removed; /* |A ∩ B| (dropped) */ + float* centroid_diff; /* c_A − c_B (dim, owned) */ + double centroid_diff_mag; + int dim; +} GeoSetDiff; +int engram_geo_setdiff(const GeoDescriptor* a, const GeoDescriptor* b, GeoSetDiff* out); +void engram_geo_setdiff_free(GeoSetDiff* s); + +/* combine(A,B): a merged descriptor — POOLED centroid + POOLED covariance + * (exact law-of-total-variance: the covariance you'd get by concatenating the two + * member clouds), re-eigendecomposed for its principal axes. Members = id-union + * (membership = max). top_axes<=0 → 8. Returns a malloc'd GeoDescriptor (free with + * engram_geo_free) in the same frame as A, or NULL on error. */ +GeoDescriptor* engram_geo_combine(const GeoDescriptor* a, const GeoDescriptor* b, + int top_axes); + +/* distance(A,B): centroid L2 + centroid cosine + closed-form Wasserstein-2 + * (Bures metric) between the two Gaussians — mirrors the proxy's _wasserstein2. */ +typedef struct { + double centroid_distance; + double centroid_cosine; + double wasserstein2; + int dim; +} GeoDistance; +int engram_geo_distance(const GeoDescriptor* a, const GeoDescriptor* b, GeoDistance* out); + +/* analogy(A,B): orthogonal PROCRUSTES transform min_R ‖A − B R‖_F, RᵀR=I (SVD) + * aligning A's principal frame to B's (extent-scaled axes, paired by rank). R is + * returned COMPACTLY as an r×r rotation within the joint axis subspace `basis` + * (r vectors of dim floats); it acts as the identity on the orthogonal complement. + * Apply it to a vector with engram_geo_analogy_apply. */ +typedef struct { + int dim; + int r; /* subspace rank; R is r×r */ + float* basis; /* r×dim row-major orthonormal basis Q (owned) */ + double* R; /* r×r rotation in Q-coords, row-major (owned) */ + double residual; /* ‖A − B R‖_F over the extent-scaled frames */ +} GeoAnalogy; +int engram_geo_analogy(const GeoDescriptor* a, const GeoDescriptor* b, GeoAnalogy* out); +/* out_vec = R·v for v ∈ R^dim: v + Σ_i (R̂c − c)_i q_i, c_i = q_i·v. dim floats. */ +void engram_geo_analogy_apply(const GeoAnalogy* an, const float* v, float* out_vec); +void engram_geo_analogy_free(GeoAnalogy* an); + /* ═══════════════════════════════════════════════════════════════════════════ * M10 — REIFICATION: densely co-wired relational neighborhoods crystallized into * FIRST-CLASS, PERSISTED store records (design doc §2; memory 885f5945). This is