will.anderson
85eee42106
Register the six engram_geo_*_json operators in the compiler builtin_arity table (bare heavy-runtime names + __ seed names, mirroring engram_activate_json) and add the engram.el module wrappers, so a compiled El (CGI) program can call them by name. The heavy-runtime C functions already existed (el_runtime.c:12287+, declared el_runtime.h:627-632); this completes the EL call surface. The shipped elc already emits a direct C call for these builtins (unknown ident-calls pass through), so no self-host compiler fold — the memory-heavy, drift-prone step — was required. Demonstrated end-to-end: a compiled geo_ops_demo.el booted a copy of the store (13,036 nodes) and produced real subtract/distance JSON on two real neighborhoods; test_geo_ops.c stays 20/20, ASan/UBSan clean. Also brace the centroid_unit normalization if/else in engram_geometry.c to clear the misleading-indentation warning (behavior-neutral).
1504 lines
69 KiB
C
1504 lines
69 KiB
C
/* engram_geometry.c — M9 FOUNDATION: relational-neighborhood geometry descriptor.
|
||
* See engram_geometry.h. Pure C11, stdlib + libm. READ-ONLY over store + vindex.
|
||
*/
|
||
#include "engram_geometry.h"
|
||
|
||
#include <stdlib.h>
|
||
#include <string.h>
|
||
#include <math.h>
|
||
#include <stdio.h>
|
||
#include <stdarg.h>
|
||
#include <time.h>
|
||
|
||
/* Must match ENGRAM_HEBB_GAIN in el_runtime.c (eff = weight*(1+GAIN*hebb)). */
|
||
#define GEO_HEBB_GAIN 0.5
|
||
/* Internal cap on the m×m Jacobi eigensolve: above this we still give centroid +
|
||
* radius but skip principal axes (honest degradation, not a lie). */
|
||
#define GEO_EIG_CAP 512
|
||
|
||
/* ───────────────────────── small dynamic member table ──────────────────────
|
||
* Neighborhoods are small (tens..few hundred), so linear-scan dedup is fine. */
|
||
typedef struct {
|
||
char** id; /* strdup'd ids */
|
||
double* memb; /* provisional membership */
|
||
float** emb; /* L2-normalized emb copy (dim floats) or NULL */
|
||
double* sal; /* stored salience */
|
||
int n, cap, dim;
|
||
} MemSet;
|
||
|
||
static int ms_init(MemSet* s, int dim){
|
||
s->n=0; s->cap=16; s->dim=dim;
|
||
s->id=calloc(s->cap,sizeof*s->id); s->memb=calloc(s->cap,sizeof*s->memb);
|
||
s->emb=calloc(s->cap,sizeof*s->emb); s->sal=calloc(s->cap,sizeof*s->sal);
|
||
return (s->id&&s->memb&&s->emb&&s->sal)?0:-1;
|
||
}
|
||
static int ms_find(const MemSet* s, const char* id){
|
||
for(int i=0;i<s->n;i++) if(strcmp(s->id[i],id)==0) return i;
|
||
return -1;
|
||
}
|
||
/* Insert or bump membership (keep the max). Returns member index or <0 on OOM. */
|
||
static int ms_upsert(MemSet* s, const char* id, double memb){
|
||
int i=ms_find(s,id);
|
||
if(i>=0){ if(memb>s->memb[i]) s->memb[i]=memb; return i; }
|
||
if(s->n==s->cap){
|
||
int nc=s->cap*2;
|
||
char** a=realloc(s->id,nc*sizeof*a); if(!a) return -1; s->id=a;
|
||
double* b=realloc(s->memb,nc*sizeof*b); if(!b) return -1; s->memb=b;
|
||
float** c=realloc(s->emb,nc*sizeof*c); if(!c) return -1; s->emb=c;
|
||
double* d=realloc(s->sal,nc*sizeof*d); if(!d) return -1; s->sal=d;
|
||
s->cap=nc;
|
||
}
|
||
s->id[s->n]=strdup(id); if(!s->id[s->n]) return -1;
|
||
s->memb[s->n]=memb; s->emb[s->n]=NULL; s->sal[s->n]=0.0;
|
||
return s->n++;
|
||
}
|
||
static void ms_free(MemSet* s){
|
||
for(int i=0;i<s->n;i++){ free(s->id[i]); free(s->emb[i]); }
|
||
free(s->id); free(s->memb); free(s->emb); free(s->sal);
|
||
}
|
||
|
||
/* L2-normalize a copy of v into out (dim floats). Returns 0, or -1 if ~zero. */
|
||
static int normcopy(const float* v, int dim, float* out){
|
||
double s=0; for(int i=0;i<dim;i++) s+=(double)v[i]*v[i];
|
||
double n=sqrt(s); if(n<1e-12) return -1;
|
||
for(int i=0;i<dim;i++) out[i]=(float)(v[i]/n);
|
||
return 0;
|
||
}
|
||
/* ── mean-centered cosine geometry (whitening the anisotropic emb space) ──────
|
||
* All three operate on RAW (unnormalized-here, but member embs are already unit)
|
||
* vectors, subtracting the global-mean offset `gm` on the fly. With gm all-zeros
|
||
* they reduce EXACTLY to the raw unit-space cosine — so the same code path serves
|
||
* both the centered and legacy-raw modes. */
|
||
static double cnorm2(const float* a, const float* gm, int dim){
|
||
double s=0; for(int d=0;d<dim;d++){ double v=(double)a[d]-gm[d]; s+=v*v; } return s;
|
||
}
|
||
static double cdot(const float* a, const float* b, const float* gm, int dim){
|
||
double s=0; for(int d=0;d<dim;d++){ double av=(double)a[d]-gm[d], bv=(double)b[d]-gm[d]; s+=av*bv; }
|
||
return s;
|
||
}
|
||
static double ccos(const float* a, const float* b, const float* gm, int dim){
|
||
double na=sqrt(cnorm2(a,gm,dim)), nb=sqrt(cnorm2(b,gm,dim));
|
||
if(na<1e-12||nb<1e-12) return 0.0;
|
||
double c=cdot(a,b,gm,dim)/(na*nb); if(c>1)c=1; if(c<-1)c=-1; return c;
|
||
}
|
||
/* cosine of (a-gm) against a pre-centered UNIT direction `dir`. */
|
||
static double ccos_dir(const float* a, const float* gm, const float* dir, int dim){
|
||
double na=sqrt(cnorm2(a,gm,dim)); if(na<1e-12) return 0.0;
|
||
double s=0; for(int d=0;d<dim;d++) s+=((double)a[d]-gm[d])*(double)dir[d];
|
||
double c=s/na; if(c>1)c=1; if(c<-1)c=-1; return c;
|
||
}
|
||
|
||
/* ───────────────────────── global-mean cache ───────────────────────────────
|
||
* Store-derived centering offset: the mean of the L2-normalized embeddings over
|
||
* the embed-eligible set. See engram_geometry.h for the anisotropy rationale. */
|
||
struct GeoMeanCache { float* mean; int dim; uint64_t n; };
|
||
|
||
typedef struct { double* sum; int dim; uint64_t n; int err; } GeoMeanAcc;
|
||
/* The reified records (Neighborhood / GeoMeanFrame) carry an emb (centroid / mean)
|
||
* but are STRUCTURE, not corpus content — they must never pollute the store-wide
|
||
* mean, the hub scan, or the descriptor. One predicate, used everywhere. */
|
||
static int geo_is_reified_type(const char* nt){
|
||
return nt && (strcmp(nt,ENGRAM_GEO_NBHD_TYPE)==0 ||
|
||
strcmp(nt,ENGRAM_GEO_MEANFRAME_TYPE)==0);
|
||
}
|
||
/* Identify a structural record by its id convention (no store read needed), so the
|
||
* descriptor never admits a reified Neighborhood / GeoMeanFrame as a neighborhood
|
||
* MEMBER even when the ANN index or adjacency still references it (re-reify/refresh
|
||
* on a store that already holds reified records; ad-hoc descriptors alike). */
|
||
static int geo_is_structural_id(const char* id){
|
||
if(!id) return 0;
|
||
if(strcmp(id,ENGRAM_GEO_MEANFRAME_ID)==0) return 1;
|
||
size_t p=strlen(ENGRAM_GEO_NBHD_ID_PREFIX);
|
||
return strncmp(id,ENGRAM_GEO_NBHD_ID_PREFIX,p)==0;
|
||
}
|
||
static void geo_mean_cb(const StoreNode* n, void* ctx){
|
||
GeoMeanAcc* a=ctx; if(a->err) return;
|
||
if(geo_is_reified_type(n->node_type)) return; /* skip structural records */
|
||
if(!(n->emb && n->emb_dim>0)) return; /* skip unembedded */
|
||
if(a->dim==0){
|
||
a->dim=n->emb_dim;
|
||
a->sum=calloc((size_t)a->dim,sizeof(double));
|
||
if(!a->sum){ a->err=1; return; }
|
||
}
|
||
if(n->emb_dim!=a->dim) return; /* skip off-dim */
|
||
double s=0; for(int d=0;d<a->dim;d++) s+=(double)n->emb[d]*n->emb[d];
|
||
double nn=sqrt(s); if(nn<1e-12) return; /* skip ~zero */
|
||
for(int d=0;d<a->dim;d++) a->sum[d]+=(double)n->emb[d]/nn;
|
||
a->n++;
|
||
}
|
||
typedef struct { uint64_t n; } GeoCntAcc;
|
||
static void geo_cnt_cb(const StoreNode* n, void* ctx){
|
||
if(n->emb && n->emb_dim>0) ((GeoCntAcc*)ctx)->n++;
|
||
}
|
||
|
||
GeoMeanCache* engram_geo_mean_build(EngramPagedStore* store){
|
||
if(!store) return NULL;
|
||
GeoMeanAcc a; memset(&a,0,sizeof a);
|
||
if(store_scan_nodes(store,geo_mean_cb,&a)<0){ free(a.sum); return NULL; }
|
||
if(a.err || a.n==0 || !a.sum){ free(a.sum); return NULL; }
|
||
GeoMeanCache* c=calloc(1,sizeof*c);
|
||
if(!c){ free(a.sum); return NULL; }
|
||
c->mean=malloc((size_t)a.dim*sizeof(float));
|
||
if(!c->mean){ free(a.sum); free(c); return NULL; }
|
||
for(int d=0;d<a.dim;d++) c->mean[d]=(float)(a.sum[d]/(double)a.n);
|
||
c->dim=a.dim; c->n=a.n; free(a.sum);
|
||
return c;
|
||
}
|
||
const float* engram_geo_mean_vec(const GeoMeanCache* c){ return c?c->mean:NULL; }
|
||
int engram_geo_mean_dim(const GeoMeanCache* c){ return c?c->dim:0; }
|
||
uint64_t engram_geo_mean_count(const GeoMeanCache* c){ return c?c->n:0; }
|
||
|
||
int engram_geo_mean_maybe_refresh(GeoMeanCache* c, EngramPagedStore* store, double frac){
|
||
if(!c||!store) return -1;
|
||
GeoCntAcc cn={0};
|
||
if(store_scan_nodes(store,geo_cnt_cb,&cn)<0) return -1;
|
||
double base=(double)(c->n?c->n:1);
|
||
double drift=fabs((double)cn.n-(double)c->n)/base;
|
||
if(drift<=frac) return 0; /* no significant change */
|
||
GeoMeanAcc a; memset(&a,0,sizeof a);
|
||
if(store_scan_nodes(store,geo_mean_cb,&a)<0){ free(a.sum); return -1; }
|
||
if(a.err || a.n==0 || !a.sum){ free(a.sum); return -1; }
|
||
float* nm=malloc((size_t)a.dim*sizeof(float));
|
||
if(!nm){ free(a.sum); return -1; }
|
||
for(int d=0;d<a.dim;d++) nm[d]=(float)(a.sum[d]/(double)a.n);
|
||
free(a.sum); free(c->mean);
|
||
c->mean=nm; c->dim=a.dim; c->n=a.n;
|
||
return 1;
|
||
}
|
||
void engram_geo_mean_free(GeoMeanCache* c){ if(c){ free(c->mean); free(c); } }
|
||
|
||
/* Attach a member's emb (normalized) + salience by point-reading the store. */
|
||
static void ms_load_node(MemSet* s, int i, EngramPagedStore* st){
|
||
StoreNode nn; memset(&nn,0,sizeof nn);
|
||
if(store_get_node(st, s->id[i], &nn)!=1){ return; }
|
||
s->sal[i]=nn.salience;
|
||
if(nn.emb && nn.emb_dim==s->dim){
|
||
float* e=malloc((size_t)s->dim*sizeof(float));
|
||
if(e && normcopy(nn.emb,s->dim,e)==0) s->emb[i]=e; else free(e);
|
||
}
|
||
store_node_free(&nn);
|
||
}
|
||
|
||
/* ───────────────────────── Jacobi symmetric eigensolver ─────────────────────
|
||
* Cyclic Jacobi on a dense symmetric m×m matrix A (row-major, overwritten).
|
||
* Eigenvalues -> w[m]; eigenvectors (columns) -> V[m*m]. Robust, libm-only. */
|
||
static void jacobi_sym(double* A, int m, double* w, double* V){
|
||
for(int i=0;i<m;i++){ for(int j=0;j<m;j++) V[i*m+j]=(i==j)?1.0:0.0; }
|
||
for(int sweep=0; sweep<100; sweep++){
|
||
double off=0; for(int p=0;p<m;p++) for(int q=p+1;q<m;q++) off+=A[p*m+q]*A[p*m+q];
|
||
if(off < 1e-18) break;
|
||
for(int p=0;p<m;p++) for(int q=p+1;q<m;q++){
|
||
double apq=A[p*m+q]; if(fabs(apq)<1e-300) continue;
|
||
double app=A[p*m+p], aqq=A[q*m+q];
|
||
double phi=0.5*atan2(2*apq, aqq-app);
|
||
double c=cos(phi), sn=sin(phi);
|
||
for(int k=0;k<m;k++){
|
||
double akp=A[k*m+p], akq=A[k*m+q];
|
||
A[k*m+p]=c*akp - sn*akq; A[k*m+q]=sn*akp + c*akq;
|
||
}
|
||
for(int k=0;k<m;k++){
|
||
double apk=A[p*m+k], aqk=A[q*m+k];
|
||
A[p*m+k]=c*apk - sn*aqk; A[q*m+k]=sn*apk + c*aqk;
|
||
}
|
||
for(int k=0;k<m;k++){
|
||
double vkp=V[k*m+p], vkq=V[k*m+q];
|
||
V[k*m+p]=c*vkp - sn*vkq; V[k*m+q]=sn*vkp + c*vkq;
|
||
}
|
||
}
|
||
}
|
||
for(int i=0;i<m;i++) w[i]=A[i*m+i];
|
||
}
|
||
|
||
void engram_geo_default_params(GeoParams* p){
|
||
if(!p) return;
|
||
p->ann_k=24; p->hop_relational=1; p->edge_min_weight=0.05;
|
||
p->kcore_k=0; p->top_axes=8; p->max_members=400;
|
||
}
|
||
|
||
/* Effective hebb-weighted edge strength, matching eg_edge_eff_weight. */
|
||
static double eff_w(double weight, double hebb){
|
||
double w = weight * (1.0 + GEO_HEBB_GAIN*hebb);
|
||
if(w>1.0) w=1.0; if(w<0.0) w=0.0; return w;
|
||
}
|
||
|
||
GeoDescriptor* engram_geometry_descriptor(
|
||
EngramPagedStore* store, VIndex* vindex,
|
||
char** vids, int n_vids,
|
||
const char* const* seed_ids, size_t n_seeds,
|
||
const GeoParams* params,
|
||
const float* global_mean)
|
||
{
|
||
if(!store || !seed_ids || n_seeds==0) return NULL;
|
||
GeoParams P; if(params) P=*params; else engram_geo_default_params(&P);
|
||
|
||
int dim = 0;
|
||
/* infer dim from the first embedded seed */
|
||
for(size_t i=0;i<n_seeds && dim==0;i++){
|
||
StoreNode nn; memset(&nn,0,sizeof nn);
|
||
if(store_get_node(store, seed_ids[i], &nn)==1){
|
||
if(nn.emb && nn.emb_dim>0) dim=nn.emb_dim;
|
||
}
|
||
store_node_free(&nn);
|
||
}
|
||
if(dim==0) dim = 768; /* no embedded seed: still build the relational side */
|
||
|
||
/* Centering offset. When a global_mean is supplied the semantic cosine math
|
||
* runs in mean-centered (isotropic) space; otherwise GM is an all-zeros
|
||
* vector so the identical code path reproduces raw unit-space cosines. */
|
||
int centered = (global_mean != NULL);
|
||
float* zeros = NULL;
|
||
const float* GM;
|
||
if(centered) GM = global_mean;
|
||
else { zeros = calloc((size_t)dim,sizeof(float));
|
||
if(!zeros) return NULL; GM = zeros; }
|
||
|
||
MemSet ms; if(ms_init(&ms,dim)!=0){ ms_free(&ms); free(zeros); return NULL; }
|
||
|
||
/* 1. seeds (membership 1.0) */
|
||
for(size_t i=0;i<n_seeds;i++) ms_upsert(&ms, seed_ids[i], 1.0);
|
||
int n_seed_members = ms.n;
|
||
for(int i=0;i<ms.n;i++) ms_load_node(&ms,i,store);
|
||
|
||
/* provisional centroid from seed embeddings (for the ANN query) */
|
||
float* prov = calloc((size_t)dim,sizeof(float));
|
||
int prov_n=0;
|
||
for(int i=0;i<n_seed_members;i++) if(ms.emb[i]){
|
||
for(int d=0;d<dim;d++) prov[d]+=ms.emb[i][d]; prov_n++;
|
||
}
|
||
if(prov_n){ for(int d=0;d<dim;d++) prov[d]/=(float)prov_n; }
|
||
|
||
/* 2. semantic expansion via vindex ANN around the provisional centroid */
|
||
if(vindex && vids && P.ann_k>0 && prov_n>0){
|
||
int k=P.ann_k*(int)n_seeds; if(k<P.ann_k) k=P.ann_k; if(k>n_vids) k=n_vids;
|
||
uint64_t* rids=malloc((size_t)k*sizeof(uint64_t));
|
||
float* dd=malloc((size_t)k*sizeof(float));
|
||
if(rids&&dd){
|
||
int got=vindex_search(vindex, prov, k, 0, rids, dd);
|
||
for(int r=0;r<got;r++){
|
||
if(rids[r]>=(uint64_t)n_vids) continue;
|
||
if(geo_is_structural_id(vids[rids[r]])) continue; /* never a member */
|
||
double memb = 1.0 - (double)dd[r]; /* cosine sim in [-1,1] */
|
||
if(memb<0) memb=0;
|
||
int mi=ms_upsert(&ms, vids[rids[r]], memb*0.9); /* <1: not a seed */
|
||
if(mi>=0 && !ms.emb[mi]) ms_load_node(&ms,mi,store);
|
||
}
|
||
}
|
||
free(rids); free(dd);
|
||
}
|
||
free(prov);
|
||
|
||
/* 3. relational expansion: seeds' hebb neighbors become members */
|
||
if(P.hop_relational){
|
||
for(int i=0;i<n_seed_members;i++){
|
||
StoreEdge* es=NULL; size_t ne=0;
|
||
if(store_get_edges_from(store, ms.id[i], &es, &ne)==0 && es){
|
||
for(size_t e=0;e<ne;e++){
|
||
if(es[e].tombstoned || es[e].inhibitory) continue;
|
||
if(es[e].relation && strcmp(es[e].relation,ENGRAM_GEO_MEMBER_RELATION)==0) continue;
|
||
if(geo_is_structural_id(es[e].to_id)) continue;
|
||
double w=eff_w(es[e].weight, es[e].hebb);
|
||
if(w < P.edge_min_weight) continue;
|
||
int mi=ms_upsert(&ms, es[e].to_id, w);
|
||
if(mi>=0 && !ms.emb[mi]) ms_load_node(&ms,mi,store);
|
||
}
|
||
}
|
||
store_edges_free(es,ne);
|
||
es=NULL; ne=0;
|
||
if(store_get_edges_to(store, ms.id[i], &es, &ne)==0 && es){
|
||
for(size_t e=0;e<ne;e++){
|
||
if(es[e].tombstoned || es[e].inhibitory) continue;
|
||
if(es[e].relation && strcmp(es[e].relation,ENGRAM_GEO_MEMBER_RELATION)==0) continue;
|
||
if(geo_is_structural_id(es[e].from_id)) continue;
|
||
double w=eff_w(es[e].weight, es[e].hebb);
|
||
if(w < P.edge_min_weight) continue;
|
||
int mi=ms_upsert(&ms, es[e].from_id, w);
|
||
if(mi>=0 && !ms.emb[mi]) ms_load_node(&ms,mi,store);
|
||
}
|
||
}
|
||
store_edges_free(es,ne);
|
||
}
|
||
}
|
||
|
||
/* optional cap: keep the highest-membership members (guards eigensolve) */
|
||
if(P.max_members>0 && ms.n>P.max_members){
|
||
/* simple selection: repeatedly drop the min-membership non-seed member */
|
||
while(ms.n>P.max_members){
|
||
int worst=-1; double wv=1e30;
|
||
for(int i=n_seed_members;i<ms.n;i++) if(ms.memb[i]<wv){wv=ms.memb[i];worst=i;}
|
||
if(worst<0) break;
|
||
free(ms.id[worst]); free(ms.emb[worst]);
|
||
ms.id[worst]=ms.id[ms.n-1]; ms.emb[worst]=ms.emb[ms.n-1];
|
||
ms.memb[worst]=ms.memb[ms.n-1]; ms.sal[worst]=ms.sal[ms.n-1];
|
||
ms.n--;
|
||
}
|
||
}
|
||
|
||
int M = ms.n;
|
||
/* ── final centroid over all embedded members ── */
|
||
float* centroid=calloc((size_t)dim,sizeof(float));
|
||
int nemb=0; int* eidx=malloc((size_t)M*sizeof(int));
|
||
for(int i=0;i<M;i++) if(ms.emb[i]){ eidx[nemb++]=i;
|
||
for(int d=0;d<dim;d++) centroid[d]+=ms.emb[i][d]; }
|
||
if(nemb){ for(int d=0;d<dim;d++) centroid[d]/=(float)nemb; }
|
||
|
||
/* ── radius + per-member cosine distance to centroid (CENTERED frame) ── */
|
||
double total_var=0;
|
||
double* distc=calloc((size_t)M,sizeof(double));
|
||
/* centered centroid (= raw centroid - global mean) and its unit direction */
|
||
float* ccen=malloc((size_t)dim*sizeof(float));
|
||
for(int d=0;d<dim;d++) ccen[d]=centroid[d]-GM[d];
|
||
float* cdir=malloc((size_t)dim*sizeof(float));
|
||
int have_cdir = (nemb>0 && normcopy(ccen,dim,cdir)==0);
|
||
for(int i=0;i<M;i++){
|
||
if(ms.emb[i] && have_cdir){
|
||
double cs=ccos_dir(ms.emb[i],GM,cdir,dim);
|
||
distc[i]=1.0-cs;
|
||
} else distc[i]=-1.0; /* unknown */
|
||
}
|
||
/* variance = mean squared Euclid distance of normalized embs to centroid */
|
||
for(int j=0;j<nemb;j++){
|
||
int i=eidx[j]; double s=0;
|
||
for(int d=0;d<dim;d++){ double df=(double)ms.emb[i][d]-centroid[d]; s+=df*df; }
|
||
total_var+=s;
|
||
}
|
||
if(nemb) total_var/=nemb;
|
||
double radius=sqrt(total_var>0?total_var:0);
|
||
|
||
/* ── principal axes via dual PCA (Jacobi on the m×m Gram of centered embs) ──
|
||
* Skipped entirely when top_axes==0: the eigensolve is the dominant cost, and
|
||
* priming needs only members+membership, so reified records that don't want the
|
||
* ellipsoid pass top_axes=0 and pay nothing here (centroid+radius still filled). */
|
||
int n_axes=0; GeoAxis* axes=NULL;
|
||
if(P.top_axes>0 && nemb>=2 && nemb<=GEO_EIG_CAP){
|
||
int m=nemb;
|
||
/* centered, row-major m×dim */
|
||
float* Xc=malloc((size_t)m*dim*sizeof(float));
|
||
for(int j=0;j<m;j++){ int i=eidx[j];
|
||
for(int d=0;d<dim;d++) Xc[(size_t)j*dim+d]=ms.emb[i][d]-centroid[d]; }
|
||
double* G=malloc((size_t)m*m*sizeof(double));
|
||
for(int a=0;a<m;a++) for(int b=a;b<m;b++){
|
||
double s=0; for(int d=0;d<dim;d++) s+=(double)Xc[(size_t)a*dim+d]*Xc[(size_t)b*dim+d];
|
||
G[a*m+b]=s; G[b*m+a]=s;
|
||
}
|
||
double* w=malloc((size_t)m*sizeof(double));
|
||
double* V=malloc((size_t)m*m*sizeof(double));
|
||
jacobi_sym(G,m,w,V);
|
||
/* sort eigenvalue indices descending */
|
||
int* ord=malloc((size_t)m*sizeof(int));
|
||
for(int i=0;i<m;i++) ord[i]=i;
|
||
for(int a=0;a<m;a++) for(int b=a+1;b<m;b++) if(w[ord[b]]>w[ord[a]]){int t=ord[a];ord[a]=ord[b];ord[b]=t;}
|
||
int keep=P.top_axes; if(keep>m) keep=m; if(keep<0) keep=0;
|
||
axes=calloc((size_t)keep,sizeof(GeoAxis));
|
||
for(int t=0;t<keep;t++){
|
||
int c=ord[t];
|
||
double lam=w[c]; if(lam<0) lam=0;
|
||
double eigcov = lam/(double)(m-1); /* covariance eigenvalue */
|
||
/* principal axis in R^dim: a = Xc^T u_c, then unit-normalize */
|
||
float* ax=calloc((size_t)dim,sizeof(float));
|
||
for(int d=0;d<dim;d++){ double s=0;
|
||
for(int j=0;j<m;j++) s+=(double)V[j*m+c]*Xc[(size_t)j*dim+d];
|
||
ax[d]=(float)s; }
|
||
double nn=0; for(int d=0;d<dim;d++) nn+=(double)ax[d]*ax[d]; nn=sqrt(nn);
|
||
if(nn>1e-12) for(int d=0;d<dim;d++) ax[d]=(float)(ax[d]/nn);
|
||
axes[t].axis=ax; axes[t].extent=sqrt(eigcov);
|
||
n_axes++;
|
||
}
|
||
free(Xc); free(G); free(w); free(V); free(ord);
|
||
}
|
||
|
||
/* ── skeleton: internal hebb edges among members + centrality gradient ── */
|
||
GeoEdge* edges=NULL; int n_edges=0, cap_e=0;
|
||
double* centrality=calloc((size_t)M,sizeof(double));
|
||
int* degree=calloc((size_t)M,sizeof(int));
|
||
/* co-registration accumulators */
|
||
double cr_n=0, cr_sx=0, cr_sy=0, cr_sxx=0, cr_syy=0, cr_sxy=0;
|
||
for(int i=0;i<M;i++){
|
||
StoreEdge* es=NULL; size_t ne=0;
|
||
if(store_get_edges_from(store, ms.id[i], &es, &ne)==0 && es){
|
||
for(size_t e=0;e<ne;e++){
|
||
if(es[e].tombstoned || es[e].inhibitory) continue;
|
||
int j=ms_find(&ms, es[e].to_id);
|
||
if(j<0 || j<=i) continue; /* internal, undirected, i<j only */
|
||
double w=eff_w(es[e].weight, es[e].hebb);
|
||
if(w < P.edge_min_weight) continue;
|
||
if(n_edges==cap_e){ cap_e=cap_e?cap_e*2:32;
|
||
GeoEdge* t=realloc(edges,(size_t)cap_e*sizeof(GeoEdge)); if(!t) break; edges=t; }
|
||
edges[n_edges].a=(uint32_t)i; edges[n_edges].b=(uint32_t)j;
|
||
edges[n_edges].eff_weight=w; edges[n_edges].hebb=es[e].hebb;
|
||
n_edges++;
|
||
centrality[i]+=w; centrality[j]+=w; degree[i]++; degree[j]++;
|
||
/* co-registration: relational strength vs semantic proximity
|
||
* (semantic proximity measured in the CENTERED frame). */
|
||
if(ms.emb[i] && ms.emb[j]){
|
||
double cs=ccos(ms.emb[i],ms.emb[j],GM,dim);
|
||
double x=w, y=cs;
|
||
cr_n++; cr_sx+=x; cr_sy+=y; cr_sxx+=x*x; cr_syy+=y*y; cr_sxy+=x*y;
|
||
}
|
||
}
|
||
}
|
||
store_edges_free(es,ne);
|
||
}
|
||
double co_reg=0;
|
||
if(cr_n>=2){
|
||
double cov=cr_sxy - cr_sx*cr_sy/cr_n;
|
||
double vx=cr_sxx - cr_sx*cr_sx/cr_n, vy=cr_syy - cr_sy*cr_sy/cr_n;
|
||
if(vx>1e-12 && vy>1e-12) co_reg=cov/sqrt(vx*vy);
|
||
}
|
||
|
||
/* ── k-core: peel members by internal degree to get core numbers ── */
|
||
int* core=calloc((size_t)M,sizeof(int));
|
||
{
|
||
int* deg=malloc((size_t)M*sizeof(int));
|
||
int* removed=calloc((size_t)M,sizeof(int));
|
||
for(int i=0;i<M;i++) deg[i]=degree[i];
|
||
int level=0, remaining=M;
|
||
while(remaining>0){
|
||
int progressed=0;
|
||
for(int i=0;i<M;i++){
|
||
if(!removed[i] && deg[i]<=level){
|
||
core[i]=level; removed[i]=1; remaining--; progressed=1;
|
||
/* decrement neighbors' working degree */
|
||
for(int e=0;e<n_edges;e++){
|
||
int o=-1;
|
||
if((int)edges[e].a==i && !removed[edges[e].b]) o=edges[e].b;
|
||
else if((int)edges[e].b==i && !removed[edges[e].a]) o=edges[e].a;
|
||
if(o>=0) deg[o]--;
|
||
}
|
||
}
|
||
}
|
||
if(!progressed) level++;
|
||
}
|
||
free(deg); free(removed);
|
||
}
|
||
int k_core=0; for(int i=0;i<M;i++) if(core[i]>k_core) k_core=core[i];
|
||
|
||
/* hub = highest centrality (tie-break salience) */
|
||
int hub=-1; double hv=-1;
|
||
for(int i=0;i<M;i++){ double v=centrality[i]+1e-6*ms.sal[i];
|
||
if(v>hv){hv=v;hub=i;} }
|
||
if(hub<0) hub=0;
|
||
|
||
/* ── assemble descriptor ── */
|
||
GeoDescriptor* g=calloc(1,sizeof(GeoDescriptor));
|
||
g->dim=dim;
|
||
g->hub_id = strdup(ms.id[hub]);
|
||
g->centroid = ccen; /* CENTERED centroid; transfer ownership */
|
||
free(centroid);
|
||
if(centered){
|
||
g->global_mean = malloc((size_t)dim*sizeof(float));
|
||
if(g->global_mean) memcpy(g->global_mean, GM, (size_t)dim*sizeof(float));
|
||
} else g->global_mean = NULL;
|
||
g->n_axes=n_axes; g->axes=axes;
|
||
g->total_variance=total_var; g->radius=radius;
|
||
g->n_members=M; g->n_embedded=nemb;
|
||
g->members=calloc((size_t)M,sizeof(GeoMember));
|
||
for(int i=0;i<M;i++){
|
||
g->members[i].id=strdup(ms.id[i]);
|
||
g->members[i].membership=ms.memb[i];
|
||
g->members[i].centrality=centrality[i];
|
||
g->members[i].salience=ms.sal[i];
|
||
g->members[i].core=core[i];
|
||
g->members[i].dist_centroid=distc[i];
|
||
g->members[i].embedded=ms.emb[i]?1:0;
|
||
}
|
||
g->n_edges=n_edges; g->edges=edges;
|
||
g->k_core=(P.kcore_k>0?P.kcore_k:k_core);
|
||
g->co_registration=co_reg;
|
||
|
||
free(centrality); free(degree); free(core); free(distc); free(eidx);
|
||
free(cdir); free(zeros);
|
||
ms_free(&ms);
|
||
return g;
|
||
}
|
||
|
||
void engram_geo_free(GeoDescriptor* g){
|
||
if(!g) return;
|
||
free(g->hub_id); free(g->centroid); free(g->global_mean);
|
||
for(int i=0;i<g->n_axes;i++) free(g->axes[i].axis);
|
||
free(g->axes);
|
||
for(int i=0;i<g->n_members;i++) free(g->members[i].id);
|
||
free(g->members); free(g->edges);
|
||
free(g);
|
||
}
|
||
|
||
/* ═══════════════════════════════════════════════════════════════════════════
|
||
* M-INTEROCEPTION P3 — DRIFT SENSOR primitive (descriptor displacement).
|
||
* Read-only. Measures how far descriptor B has drifted from a baseline A and
|
||
* decomposes it into GROWTH (periphery extends, core fixed) vs CORRUPTION (the
|
||
* invariant core displaces). The core is the top `core_frac` of A's members by
|
||
* centrality; the periphery is the rest. Per shared member (matched by id), the
|
||
* displacement is the change in its radial position (dist_centroid) between A
|
||
* and B; centroid separation + radius delta give the aggregate move.
|
||
*
|
||
* PREREQUISITE FLAGGED (honesty rail, design §2/§6): a LIVE self-drift reading
|
||
* needs a persisted SelfAnchor baseline descriptor to compare "now" against.
|
||
* That anchor does NOT exist yet — there is no persisted self node / anchored
|
||
* self-neighborhood in this store. This primitive therefore takes an EXPLICIT
|
||
* baseline so it is real and testable today; capturing a durable SelfAnchor
|
||
* snapshot and wiring the ENGRAM_DRIFT_SENSOR live reading is a follow-up. We do
|
||
* NOT fabricate a self silently.
|
||
* ═══════════════════════════════════════════════════════════════════════════ */
|
||
static double eg_geo_l2(const float* x, const float* y, int dim){
|
||
double s=0; for(int i=0;i<dim;i++){ double d=(double)x[i]-(double)y[i]; s+=d*d; } return sqrt(s);
|
||
}
|
||
static double eg_geo_cosv(const float* x, const float* y, int dim){
|
||
double dot=0,nx=0,ny=0;
|
||
for(int i=0;i<dim;i++){ dot+=(double)x[i]*y[i]; nx+=(double)x[i]*x[i]; ny+=(double)y[i]*y[i]; }
|
||
if(nx<=0.0||ny<=0.0) return 0.0;
|
||
return dot/(sqrt(nx)*sqrt(ny));
|
||
}
|
||
void engram_geo_displacement(const GeoDescriptor* a, const GeoDescriptor* b,
|
||
double core_frac, GeoDisplacement* out){
|
||
if(!out) return;
|
||
memset(out,0,sizeof(*out));
|
||
if(!a||!b) return;
|
||
if(a->centroid && b->centroid && a->dim==b->dim && a->dim>0){
|
||
out->centroid_sep = eg_geo_l2(a->centroid,b->centroid,a->dim);
|
||
out->centroid_cos = 1.0 - eg_geo_cosv(a->centroid,b->centroid,a->dim);
|
||
}
|
||
out->radius_delta = fabs(a->radius - b->radius);
|
||
if(!(core_frac>0.0 && core_frac<=1.0)) core_frac=0.3;
|
||
int na=a->n_members;
|
||
if(na<=0) return;
|
||
int* order=malloc((size_t)na*sizeof(int));
|
||
if(!order) return;
|
||
for(int i=0;i<na;i++) order[i]=i;
|
||
/* insertion sort by centrality desc (neighborhoods are small) */
|
||
for(int i=1;i<na;i++){ int k=order[i]; int j=i-1;
|
||
while(j>=0 && a->members[order[j]].centrality < a->members[k].centrality){ order[j+1]=order[j]; j--; }
|
||
order[j+1]=k; }
|
||
int ncore=(int)(core_frac*na+0.5); if(ncore<1) ncore=1; if(ncore>na) ncore=na;
|
||
double core_sum=0, periph_sum=0; int core_n=0, periph_n=0;
|
||
for(int r=0;r<na;r++){
|
||
const GeoMember* ma=&a->members[order[r]];
|
||
const GeoMember* mb=NULL;
|
||
for(int j=0;j<b->n_members;j++){
|
||
if(b->members[j].id && ma->id && strcmp(b->members[j].id,ma->id)==0){ mb=&b->members[j]; break; }
|
||
}
|
||
if(!mb) continue;
|
||
double disp=fabs(ma->dist_centroid - mb->dist_centroid);
|
||
if(r<ncore){ core_sum+=disp; core_n++; } else { periph_sum+=disp; periph_n++; }
|
||
}
|
||
free(order);
|
||
out->core_matched=core_n; out->periph_matched=periph_n;
|
||
out->core_disp = core_n ? core_sum/core_n : 0.0;
|
||
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;d<dim;d++) s+=(double)a[d]*b[d]; return s;
|
||
}
|
||
/* coefficient of w along a principal axis: a_k · w. */
|
||
static double geo_axis_coef(const GeoAxis* ax, const float* w, int dim){
|
||
double s=0; for(int d=0;d<dim;d++) s+=(double)ax->axis[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;d<dim;d++) out[d]=0.0f;
|
||
for(int k=0;k<g->n_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;d<dim;d++) out[d]+=(float)(c*a[d]);
|
||
}
|
||
}
|
||
/* trace(Σ_g) = Σ_k extent_k² (the retained variance). */
|
||
static double geo_cov_trace(const GeoDescriptor* g){
|
||
double s=0; for(int k=0;k<g->n_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;i<nc;i++){
|
||
for(int d=0;d<dim;d++) v[d]=(double)cand[i][d];
|
||
for(int j=0;j<r;j++){
|
||
double dot=0; const float* qj=&Q[(size_t)j*dim];
|
||
for(int d=0;d<dim;d++) dot+=v[d]*qj[d];
|
||
for(int d=0;d<dim;d++) v[d]-=dot*qj[d];
|
||
}
|
||
double nrm=0; for(int d=0;d<dim;d++) nrm+=v[d]*v[d]; nrm=sqrt(nrm);
|
||
if(nrm>1e-6){ for(int d=0;d<dim;d++) Q[(size_t)r*dim+d]=(float)(v[d]/nrm); r++; }
|
||
}
|
||
free(v);
|
||
return r;
|
||
}
|
||
/* M[i][j] = q_iᵀ Σ_g q_j over an r×dim basis Q; M is r×r (symmetrized). */
|
||
static void geo_cov_in_basis(const GeoDescriptor* g, const float* Q, int r, int dim, double* M){
|
||
float* sq=malloc((size_t)dim*sizeof(float));
|
||
if(!sq){ for(int i=0;i<r*r;i++) M[i]=0; return; }
|
||
for(int j=0;j<r;j++){
|
||
geo_cov_apply(g,&Q[(size_t)j*dim],sq,dim);
|
||
for(int i=0;i<r;i++){
|
||
const float* qi=&Q[(size_t)i*dim];
|
||
double s=0; for(int d=0;d<dim;d++) s+=(double)qi[d]*sq[d];
|
||
M[i*r+j]=s;
|
||
}
|
||
}
|
||
free(sq);
|
||
for(int i=0;i<r;i++) for(int j=i+1;j<r;j++){ double m=0.5*(M[i*r+j]+M[j*r+i]); M[i*r+j]=M[j*r+i]=m; }
|
||
}
|
||
/* C = A·B for n×n row-major. */
|
||
static void geo_matmul(const double* A, const double* B, double* C, int n){
|
||
for(int i=0;i<n;i++) for(int j=0;j<n;j++){
|
||
double s=0; for(int k=0;k<n;k++) s+=A[i*n+k]*B[k*n+j]; C[i*n+j]=s; }
|
||
}
|
||
/* R = S^{1/2} for symmetric PSD n×n (eigenvalues clipped ≥0). Mirrors _sym_sqrt. */
|
||
static void geo_symsqrt(const double* S, double* R, int n){
|
||
double* A=malloc((size_t)n*n*sizeof(double));
|
||
double* w=malloc((size_t)n*sizeof(double));
|
||
double* V=malloc((size_t)n*n*sizeof(double));
|
||
if(!A||!w||!V){ free(A);free(w);free(V); for(int i=0;i<n*n;i++) R[i]=0; return; }
|
||
memcpy(A,S,(size_t)n*n*sizeof(double));
|
||
jacobi_sym(A,n,w,V); /* V columns = eigenvectors */
|
||
for(int i=0;i<n;i++) for(int j=0;j<n;j++){
|
||
double s=0;
|
||
for(int k=0;k<n;k++){ double sq=w[k]>0?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_members<b->n_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;i<a->n_members;i++){
|
||
const char* id=a->members[i].id; if(!id) continue;
|
||
for(int j=0;j<b->n_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;k<dim;k++){ double df=(double)a->centroid[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;k<dim;k++) out->intersection_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;i<o->n_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;k<mB;k++){
|
||
double coef=geo_axis_coef(&b->axes[k],a->centroid,dim);
|
||
Qc2 += coef*coef;
|
||
const float* ax=b->axes[k].axis;
|
||
for(int d=0;d<dim;d++) out->residual_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;k<mB;k++){
|
||
for(int j=0;j<a->n_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;j<a->n_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;k<mB;k++){
|
||
double coef=geo_axis_coef(&b->axes[k],a->axes[j].axis,dim);
|
||
const float* ax=b->axes[k].axis;
|
||
for(int d=0;d<dim;d++) r[d]-=(float)(coef*ax[d]);
|
||
}
|
||
double nrm=0; for(int d=0;d<dim;d++) nrm+=(double)r[d]*r[d]; nrm=sqrt(nrm);
|
||
if(nrm>1e-6){ for(int d=0;d<dim;d++) r[d]=(float)(r[d]/nrm);
|
||
out->axes[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;d<dim;d++){ float df=a->centroid[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;i<r->n_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;i<a->n_members;i++){
|
||
const char* id=a->members[i].id; if(!id) continue;
|
||
int in=0;
|
||
for(int j=0;j<b->n_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;d<dim;d++){ float df=a->centroid[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;i<s->n_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;i<dim;i++){
|
||
g->centroid[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;k<a->n_axes;k++) cand[ci++]=a->axes[k].axis;
|
||
for(int k=0;k<b->n_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;j<r;j++){
|
||
const float* qj=&Q[(size_t)j*dim];
|
||
geo_cov_apply(a,qj,sqa,dim);
|
||
geo_cov_apply(b,qj,sqb,dim);
|
||
double dq = d? geo_vdot(d,qj,dim) : 0.0;
|
||
for(int i=0;i<r;i++){
|
||
const float* qi=&Q[(size_t)i*dim];
|
||
double s=0;
|
||
for(int k=0;k<dim;k++)
|
||
s+=(double)qi[k]*(wA*sqa[k]+wB*sqb[k]);
|
||
if(d) s+=cross*dq*geo_vdot(d,qi,dim);
|
||
M[i*r+j]=s;
|
||
}
|
||
}
|
||
for(int i=0;i<r;i++) for(int j=i+1;j<r;j++){ double m=0.5*(M[i*r+j]+M[j*r+i]); M[i*r+j]=M[j*r+i]=m; }
|
||
double* w=malloc((size_t)r*sizeof(double));
|
||
double* V=malloc((size_t)r*r*sizeof(double));
|
||
if(w&&V){
|
||
jacobi_sym(M,r,w,V);
|
||
int* ord=malloc((size_t)r*sizeof(int));
|
||
for(int i=0;i<r;i++) ord[i]=i;
|
||
for(int i=0;i<r;i++) for(int j=i+1;j<r;j++) if(w[ord[j]]>w[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;t<keep && g->axes;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;dd<dim;dd++){ double s=0;
|
||
for(int i=0;i<r;i++) s+=V[i*r+c]*Q[(size_t)i*dim+dd];
|
||
ax[dd]=(float)s; }
|
||
double nn=0; for(int dd=0;dd<dim;dd++) nn+=(double)ax[dd]*ax[dd]; nn=sqrt(nn);
|
||
if(nn>1e-12) for(int dd=0;dd<dim;dd++) ax[dd]=(float)(ax[dd]/nn);
|
||
g->axes[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;i<a->n_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;j<b->n_members && g->members;j++){
|
||
const char* id=b->members[j].id; int found=-1;
|
||
for(int i=0;i<M;i++) if(g->members[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;i<M;i++) if(g->members[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;k<dim;k++){ double x=a->centroid[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;k<a->n_axes;k++) cand[ci++]=a->axes[k].axis;
|
||
for(int k=0;k<b->n_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;i<r;i++) trc += C1[i*r+i]+C2[i*r+i]-2.0*inner[i*r+i];
|
||
trace_term=trc;
|
||
}
|
||
free(C1);free(C2);free(s2);free(tmp);free(mid);free(inner);
|
||
}
|
||
free(Q); free(cand);
|
||
}
|
||
double w2=d2+trace_term; if(w2<0) w2=0;
|
||
out->wasserstein2=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_axes<b->n_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;t<a->n_axes;t++) cand[ci++]=a->axes[t].axis;
|
||
for(int t=0;t<b->n_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;c<k;c++){
|
||
double ea=a->axes[c].extent, eb=b->axes[c].extent;
|
||
for(int i=0;i<r;i++){
|
||
Ah[i*k+c]=ea*geo_axis_coef(&a->axes[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;i<r;i++) for(int j=0;j<r;j++){ double s=0;
|
||
for(int c=0;c<k;c++) s+=Ah[i*k+c]*Bh[j*k+c]; Mh[i*r+j]=s; }
|
||
/* S = MhᵀMh (r×r symmetric) → jacobi → V (cols), σ²=w */
|
||
double* S=calloc((size_t)r*r,sizeof(double));
|
||
for(int i=0;i<r;i++) for(int j=0;j<r;j++){ double s=0;
|
||
for(int l=0;l<r;l++) s+=Mh[l*r+i]*Mh[l*r+j]; S[i*r+j]=s; }
|
||
double* w=malloc((size_t)r*sizeof(double));
|
||
double* V=malloc((size_t)r*r*sizeof(double));
|
||
jacobi_sym(S,r,w,V);
|
||
/* U columns: u_c = Mh·V_c / σ_c (σ_c≈0 → u_c = V_c, identity on that dir). */
|
||
double* U=calloc((size_t)r*r,sizeof(double));
|
||
for(int c=0;c<r;c++){
|
||
double sig=w[c]>0?sqrt(w[c]):0.0;
|
||
if(sig>1e-9){
|
||
for(int i=0;i<r;i++){ double s=0;
|
||
for(int j=0;j<r;j++) s+=Mh[i*r+j]*V[j*r+c]; U[i*r+c]=s/sig; }
|
||
} else {
|
||
for(int i=0;i<r;i++) U[i*r+c]=V[i*r+c];
|
||
}
|
||
}
|
||
/* R̂ = U·Vᵀ */
|
||
out->R=calloc((size_t)r*r,sizeof(double));
|
||
for(int i=0;i<r;i++) for(int j=0;j<r;j++){ double s=0;
|
||
for(int c=0;c<r;c++) s+=U[i*r+c]*V[j*r+c]; out->R[i*r+j]=s; }
|
||
/* residual = ‖Ah − R̂·Bh‖_F */
|
||
double resid=0;
|
||
for(int c=0;c<k;c++) for(int i=0;i<r;i++){
|
||
double rb=0; for(int j=0;j<r;j++) rb+=out->R[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;d<dim;d++) out_vec[d]=v[d];
|
||
if(r<=0||!an->basis||!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;i<r;i++){ double s=0; const float* qi=&an->basis[(size_t)i*dim];
|
||
for(int d=0;d<dim;d++) s+=(double)qi[d]*v[d]; c[i]=s; }
|
||
for(int i=0;i<r;i++){ double s=0; for(int j=0;j<r;j++) s+=an->R[i*r+j]*c[j]; cp[i]=s; }
|
||
for(int i=0;i<r;i++){ double delta=cp[i]-c[i]; const float* qi=&an->basis[(size_t)i*dim];
|
||
for(int d=0;d<dim;d++) out_vec[d]+=(float)(delta*qi[d]); }
|
||
free(c); free(cp);
|
||
}
|
||
void engram_geo_analogy_free(GeoAnalogy* an){
|
||
if(!an) return;
|
||
free(an->basis); free(an->R);
|
||
memset(an,0,sizeof*an);
|
||
}
|
||
|
||
/* ═══════════════════════════════════════════════════════════════════════════
|
||
* M10 — REIFICATION: persist / load / lookup first-class neighborhood records.
|
||
* ═══════════════════════════════════════════════════════════════════════════ */
|
||
|
||
static int64_t geo_now_ms(void){
|
||
struct timespec ts;
|
||
if(clock_gettime(CLOCK_REALTIME,&ts)==0)
|
||
return (int64_t)ts.tv_sec*1000 + ts.tv_nsec/1000000;
|
||
return (int64_t)time(NULL)*1000;
|
||
}
|
||
|
||
void engram_geo_reify_default_params(GeoReifyParams* p){
|
||
if(!p) return;
|
||
p->min_weighted_degree=0;
|
||
p->max_neighborhoods=128;
|
||
p->cover_membership=0.5;
|
||
p->persist_member_edges=1;
|
||
engram_geo_default_params(&p->descriptor);
|
||
p->descriptor.top_axes=4; /* keep a small ellipsoid summary; cheap */
|
||
p->descriptor.max_members=256; /* reified neighborhoods stay compact */
|
||
}
|
||
|
||
/* ── tiny growable string builder ─────────────────────────────────────────── */
|
||
typedef struct { char* s; size_t n, cap; } SB;
|
||
static int sb_reserve(SB* b, size_t add){
|
||
if(b->n+add+1<=b->cap) return 0;
|
||
size_t nc=b->cap?b->cap:256; while(nc<b->n+add+1) nc*=2;
|
||
char* t=realloc(b->s,nc); if(!t) return -1; b->s=t; b->cap=nc; return 0;
|
||
}
|
||
static int sb_puts(SB* b, const char* s){
|
||
size_t l=strlen(s); if(sb_reserve(b,l)) return -1;
|
||
memcpy(b->s+b->n,s,l); b->n+=l; b->s[b->n]=0; return 0;
|
||
}
|
||
static int sb_fmt(SB* b, const char* fmt, ...){
|
||
char tmp[512]; va_list ap; va_start(ap,fmt);
|
||
int k=vsnprintf(tmp,sizeof tmp,fmt,ap); va_end(ap);
|
||
if(k<0) return -1; if(k>=(int)sizeof tmp) k=sizeof tmp-1;
|
||
return sb_puts(b,tmp);
|
||
}
|
||
|
||
/* Serialize a descriptor's DURABLE geometry into the GEO1 metadata schema.
|
||
* (The raw centroid is stored separately as the record's emb.) */
|
||
static char* geo_nbhd_metadata(const GeoDescriptor* g, const char* hub,
|
||
const char* meanid){
|
||
SB b={0};
|
||
if(sb_puts(&b,"GEO1\n")) { free(b.s); return NULL; }
|
||
sb_fmt(&b,"hub %s\n", hub?hub:"");
|
||
sb_fmt(&b,"mean %s\n", meanid?meanid:"");
|
||
sb_fmt(&b,"s %.9g %.9g %d %.9g %d %d\n",
|
||
g->radius, g->total_variance, g->k_core, g->co_registration,
|
||
g->n_embedded, g->n_members);
|
||
sb_puts(&b,"e");
|
||
for(int i=0;i<g->n_axes;i++) sb_fmt(&b," %.9g", g->axes[i].extent);
|
||
sb_puts(&b,"\n");
|
||
for(int i=0;i<g->n_members;i++){
|
||
sb_fmt(&b,"m %s %.9g %.9g %d\n",
|
||
g->members[i].id, g->members[i].membership,
|
||
g->members[i].centrality, g->members[i].core);
|
||
}
|
||
return b.s; /* caller frees */
|
||
}
|
||
|
||
/* ── string set (greedy-cover claimed ids) + string→id list (hub→old nbhd) ──── */
|
||
static uint64_t geo_djb2(const char* s){
|
||
uint64_t h=5381; for(;*s;s++) h=((h<<5)+h)^(unsigned char)*s; return h;
|
||
}
|
||
typedef struct SSNode { char* key; struct SSNode* next; } SSNode;
|
||
typedef struct { SSNode** b; size_t nb; } SSet;
|
||
static void ss_init(SSet* s, size_t nb){ s->nb=nb; s->b=calloc(nb,sizeof*s->b); }
|
||
static int ss_has(const SSet* s, const char* k){
|
||
if(!s->b) return 0; SSNode* n=s->b[geo_djb2(k)%s->nb];
|
||
for(;n;n=n->next) if(strcmp(n->key,k)==0) return 1; return 0;
|
||
}
|
||
static void ss_add(SSet* s, const char* k){
|
||
if(!s->b||ss_has(s,k)) return; size_t i=geo_djb2(k)%s->nb;
|
||
SSNode* n=malloc(sizeof*n); if(!n) return; n->key=strdup(k); n->next=s->b[i]; s->b[i]=n;
|
||
}
|
||
static void ss_free(SSet* s){
|
||
if(!s->b) return;
|
||
for(size_t i=0;i<s->nb;i++){ SSNode* n=s->b[i]; while(n){ SSNode* x=n->next; free(n->key); free(n); n=x; } }
|
||
free(s->b); s->b=NULL;
|
||
}
|
||
|
||
typedef struct { char** id; int n, cap; } StrVec;
|
||
static void sv_push(StrVec* v, const char* s){
|
||
if(v->n==v->cap){ v->cap=v->cap?v->cap*2:64; v->id=realloc(v->id,(size_t)v->cap*sizeof*v->id); }
|
||
v->id[v->n++]=strdup(s);
|
||
}
|
||
static void sv_free(StrVec* v){ for(int i=0;i<v->n;i++) free(v->id[i]); free(v->id); }
|
||
|
||
/* pass 1 collector: all non-structural node ids; also record existing Neighborhood
|
||
* records as (hub -> old_id) so a re-reify supersedes the prior version. */
|
||
typedef struct {
|
||
StrVec cand; /* candidate node ids (content nodes) */
|
||
StrVec old_hub, old_id; /* parallel: existing nbhd hub + its id */
|
||
} ReifyScan;
|
||
static void geo_reify_scan_cb(const StoreNode* n, void* ctx){
|
||
ReifyScan* rs=ctx; if(!n->id||!n->node_type) { if(n->id) sv_push(&rs->cand,n->id); return; }
|
||
if(strcmp(n->node_type,ENGRAM_GEO_NBHD_TYPE)==0){
|
||
/* parse hub from metadata GEO1 (line "hub <id>") for supersede lineage */
|
||
const char* md=n->metadata?n->metadata:"";
|
||
const char* p=strstr(md,"hub ");
|
||
if(p && (p==md || p[-1]=='\n')){
|
||
p+=4; const char* e=p; while(*e && *e!='\n') e++;
|
||
char* hub=strndup(p,(size_t)(e-p));
|
||
sv_push(&rs->old_hub,hub); sv_push(&rs->old_id,n->id); free(hub);
|
||
}
|
||
return; /* structural: not a candidate */
|
||
}
|
||
if(strcmp(n->node_type,ENGRAM_GEO_MEANFRAME_TYPE)==0) return;
|
||
sv_push(&rs->cand,n->id);
|
||
}
|
||
|
||
/* weighted strong-edge degree of a node (from+to), matching eff_w/threshold. */
|
||
static double geo_weighted_degree(EngramPagedStore* st, const char* id, double emin){
|
||
double deg=0; StoreEdge* es=NULL; size_t ne=0;
|
||
if(store_get_edges_from(st,id,&es,&ne)==0 && es){
|
||
for(size_t e=0;e<ne;e++){ if(es[e].tombstoned||es[e].inhibitory) continue;
|
||
double w=eff_w(es[e].weight,es[e].hebb); if(w>=emin) deg+=w; }
|
||
}
|
||
store_edges_free(es,ne); es=NULL; ne=0;
|
||
if(store_get_edges_to(st,id,&es,&ne)==0 && es){
|
||
for(size_t e=0;e<ne;e++){ if(es[e].tombstoned||es[e].inhibitory) continue;
|
||
double w=eff_w(es[e].weight,es[e].hebb); if(w>=emin) deg+=w; }
|
||
}
|
||
store_edges_free(es,ne);
|
||
return deg;
|
||
}
|
||
|
||
int engram_geo_reify_store(EngramPagedStore* store, VIndex* vindex,
|
||
char** vids, int n_vids,
|
||
const GeoReifyParams* params){
|
||
if(!store) return -1;
|
||
GeoReifyParams P; if(params) P=*params; else engram_geo_reify_default_params(&P);
|
||
|
||
/* 1. true store-wide mean → persist the GeoMeanFrame record (once). */
|
||
GeoMeanCache* mc=engram_geo_mean_build(store);
|
||
if(!mc) return -2;
|
||
int dim=engram_geo_mean_dim(mc);
|
||
const float* mean=engram_geo_mean_vec(mc);
|
||
int64_t now=geo_now_ms();
|
||
{ StoreNode mf; memset(&mf,0,sizeof mf);
|
||
mf.id=(char*)ENGRAM_GEO_MEANFRAME_ID; mf.node_type=(char*)ENGRAM_GEO_MEANFRAME_TYPE;
|
||
mf.content=(char*)"geo-mean-frame"; mf.tier=(char*)"Semantic"; mf.metadata=(char*)"{}";
|
||
mf.emb=(float*)mean; mf.emb_dim=dim; mf.created_at=now; mf.updated_at=now;
|
||
if(store_put_node(store,&mf)<0){ engram_geo_mean_free(mc); return -3; }
|
||
}
|
||
|
||
/* 2. scan: candidate ids + existing (hub→old id) for supersede. */
|
||
ReifyScan rs; memset(&rs,0,sizeof rs);
|
||
if(store_scan_nodes(store,geo_reify_scan_cb,&rs)<0){
|
||
sv_free(&rs.cand); sv_free(&rs.old_hub); sv_free(&rs.old_id);
|
||
engram_geo_mean_free(mc); return -4;
|
||
}
|
||
|
||
/* 3. weighted degree per candidate; sort desc. */
|
||
int N=rs.cand.n;
|
||
double* deg=malloc((size_t)N*sizeof(double));
|
||
int* ord=malloc((size_t)N*sizeof(int));
|
||
for(int i=0;i<N;i++){ deg[i]=geo_weighted_degree(store,rs.cand.id[i],P.descriptor.edge_min_weight); ord[i]=i; }
|
||
/* simple insertion-ish selection sort by degree desc (N a few thousand, one-time) */
|
||
for(int a=0;a<N;a++){ int best=a; for(int b=a+1;b<N;b++) if(deg[ord[b]]>deg[ord[best]]) best=b;
|
||
int t=ord[a]; ord[a]=ord[best]; ord[best]=t; }
|
||
|
||
/* 4. greedy non-redundant cover: reify each qualifying hub once. */
|
||
SSet claimed; ss_init(&claimed, (size_t)(N>16?N:16));
|
||
int persisted=0;
|
||
for(int oi=0; oi<N && persisted<P.max_neighborhoods; oi++){
|
||
int i=ord[oi]; const char* hub=rs.cand.id[i];
|
||
if(P.min_weighted_degree>0 && deg[i]<(double)P.min_weighted_degree) break; /* sorted: rest smaller */
|
||
if(ss_has(&claimed,hub)) continue;
|
||
const char* seeds[1]={hub};
|
||
GeoDescriptor* g=engram_geometry_descriptor(store,vindex,vids,n_vids,
|
||
seeds,1,&P.descriptor,mean);
|
||
if(!g || g->n_members<=0){ if(g) engram_geo_free(g); continue; }
|
||
/* claim members above cover threshold (incl. the hub itself) */
|
||
for(int m=0;m<g->n_members;m++)
|
||
if(g->members[m].membership>=P.cover_membership) ss_add(&claimed,g->members[m].id);
|
||
|
||
/* build record: id = nbhd-<hub>-<now>, emb = RAW centroid = centered+mean */
|
||
char nid[512]; snprintf(nid,sizeof nid,"%s%s-%lld",ENGRAM_GEO_NBHD_ID_PREFIX,hub,(long long)now);
|
||
float* raw=NULL;
|
||
if(g->n_embedded>0 && g->centroid && g->global_mean){
|
||
raw=malloc((size_t)dim*sizeof(float));
|
||
if(raw) for(int d=0;d<dim;d++) raw[d]=g->centroid[d]+g->global_mean[d];
|
||
}
|
||
char* md=geo_nbhd_metadata(g,hub,ENGRAM_GEO_MEANFRAME_ID);
|
||
StoreNode nn; memset(&nn,0,sizeof nn);
|
||
nn.id=nid; nn.node_type=(char*)ENGRAM_GEO_NBHD_TYPE;
|
||
nn.content=(char*)"reified-neighborhood"; nn.tier=(char*)"Semantic";
|
||
nn.metadata=md?md:(char*)"{}"; nn.emb=raw; nn.emb_dim=raw?dim:0;
|
||
nn.created_at=now; nn.updated_at=now;
|
||
int wrc=store_put_node(store,&nn);
|
||
free(raw); free(md);
|
||
if(wrc<0){ engram_geo_free(g); continue; }
|
||
|
||
/* provenance: supersede any prior neighborhood for this hub. */
|
||
for(int k=0;k<rs.old_hub.n;k++) if(strcmp(rs.old_hub.id[k],hub)==0){
|
||
store_supersede(store, rs.old_id.id[k], nid);
|
||
}
|
||
|
||
/* member links (durable, but inert to activation — runtime skips them). */
|
||
if(P.persist_member_edges){
|
||
for(int m=0;m<g->n_members;m++){
|
||
char eid[600]; snprintf(eid,sizeof eid,"%s->%s",nid,g->members[m].id);
|
||
StoreEdge se; memset(&se,0,sizeof se);
|
||
se.id=eid; se.from_id=nid; se.to_id=g->members[m].id;
|
||
se.relation=(char*)ENGRAM_GEO_MEMBER_RELATION;
|
||
se.metadata=(char*)"{}"; se.weight=g->members[m].membership;
|
||
se.confidence=1.0; se.created_at=now; se.updated_at=now;
|
||
store_put_edge(store,&se);
|
||
}
|
||
}
|
||
engram_geo_free(g);
|
||
persisted++;
|
||
}
|
||
|
||
ss_free(&claimed);
|
||
free(deg); free(ord);
|
||
sv_free(&rs.cand); sv_free(&rs.old_hub); sv_free(&rs.old_id);
|
||
engram_geo_mean_free(mc);
|
||
return persisted;
|
||
}
|
||
|
||
/* ═══════════════ resident loaded form + hot-path lookup ═════════════════════ */
|
||
|
||
typedef struct {
|
||
char* id;
|
||
char* hub_id;
|
||
int n_members;
|
||
char** member_ids;
|
||
double* member_w;
|
||
double radius, co_reg;
|
||
int k_core, n_embedded;
|
||
float* centroid_raw; /* dim floats or NULL */
|
||
float* centroid_unit; /* centered+normalized (finalize) or NULL */
|
||
int dim;
|
||
GeoNeighborhood view;
|
||
} RNbhd;
|
||
|
||
typedef struct RE { char* id; int nbhd; double w; struct RE* next; } RE;
|
||
|
||
struct GeoReifyIndex {
|
||
RNbhd* nb; int n, cap;
|
||
float* mean; int mean_dim;
|
||
RE** buckets; size_t nbuckets;
|
||
double* score; /* scratch[n], reused per lookup */
|
||
};
|
||
|
||
GeoReifyIndex* engram_geo_reify_index_new(void){
|
||
GeoReifyIndex* ix=calloc(1,sizeof*ix); return ix;
|
||
}
|
||
|
||
/* parse a GEO1 metadata blob into an RNbhd (members + scalars). */
|
||
static int geo_parse_nbhd(const char* md, RNbhd* r){
|
||
if(!md) return -1;
|
||
if(strncmp(md,"GEO1",4)!=0) return -1;
|
||
/* count member lines to size arrays */
|
||
int cap=0; for(const char* p=md; (p=strstr(p,"\nm ")); p+=3) cap++;
|
||
r->member_ids=cap?calloc((size_t)cap,sizeof(char*)):NULL;
|
||
r->member_w =cap?calloc((size_t)cap,sizeof(double)):NULL;
|
||
r->n_members=0;
|
||
const char* line=md;
|
||
while(line && *line){
|
||
const char* nl=strchr(line,'\n');
|
||
size_t len= nl? (size_t)(nl-line) : strlen(line);
|
||
char buf[600]; if(len>=sizeof buf) len=sizeof buf-1;
|
||
memcpy(buf,line,len); buf[len]=0;
|
||
if(buf[0]=='h'&&buf[1]=='u'&&buf[2]=='b'&&buf[3]==' '){
|
||
free(r->hub_id); r->hub_id=strdup(buf+4);
|
||
} else if(buf[0]=='s'&&buf[1]==' '){
|
||
int kc=0,ne=0,nm=0; double rad=0,tv=0,cr=0;
|
||
sscanf(buf+2,"%lf %lf %d %lf %d %d",&rad,&tv,&kc,&cr,&ne,&nm);
|
||
r->radius=rad; r->co_reg=cr; r->k_core=kc; r->n_embedded=ne;
|
||
} else if(buf[0]=='m'&&buf[1]==' '){
|
||
char mid[512]; double w=0,c=0; int core=0;
|
||
if(sscanf(buf+2,"%511s %lf %lf %d",mid,&w,&c,&core)>=2 && r->member_ids){
|
||
r->member_ids[r->n_members]=strdup(mid);
|
||
r->member_w[r->n_members]=w;
|
||
r->n_members++;
|
||
}
|
||
}
|
||
line = nl? nl+1 : NULL;
|
||
}
|
||
return 0;
|
||
}
|
||
|
||
int engram_geo_reify_index_add(GeoReifyIndex* ix, const StoreNode* n){
|
||
if(!ix||!n||!n->node_type) return 0;
|
||
if(strcmp(n->node_type,ENGRAM_GEO_MEANFRAME_TYPE)==0){
|
||
if(n->emb && n->emb_dim>0){
|
||
free(ix->mean);
|
||
ix->mean=malloc((size_t)n->emb_dim*sizeof(float));
|
||
if(ix->mean){ memcpy(ix->mean,n->emb,(size_t)n->emb_dim*sizeof(float)); ix->mean_dim=n->emb_dim; }
|
||
}
|
||
return 0;
|
||
}
|
||
if(strcmp(n->node_type,ENGRAM_GEO_NBHD_TYPE)!=0) return 0;
|
||
if(ix->n==ix->cap){ ix->cap=ix->cap?ix->cap*2:16;
|
||
RNbhd* t=realloc(ix->nb,(size_t)ix->cap*sizeof*t); if(!t) return -1; ix->nb=t; }
|
||
RNbhd* r=&ix->nb[ix->n]; memset(r,0,sizeof*r);
|
||
r->id=strdup(n->id?n->id:"");
|
||
if(geo_parse_nbhd(n->metadata,r)!=0){ free(r->id); return 0; } /* skip malformed */
|
||
if(n->emb && n->emb_dim>0){
|
||
r->dim=n->emb_dim;
|
||
r->centroid_raw=malloc((size_t)n->emb_dim*sizeof(float));
|
||
if(r->centroid_raw) memcpy(r->centroid_raw,n->emb,(size_t)n->emb_dim*sizeof(float));
|
||
}
|
||
ix->n++;
|
||
return 0;
|
||
}
|
||
|
||
int engram_geo_reify_index_finalize(GeoReifyIndex* ix){
|
||
if(!ix) return -1;
|
||
/* member → neighborhood hash */
|
||
size_t total=0; for(int i=0;i<ix->n;i++) total+=(size_t)ix->nb[i].n_members;
|
||
ix->nbuckets = total? (total*2+1) : 1;
|
||
ix->buckets=calloc(ix->nbuckets,sizeof(RE*));
|
||
if(!ix->buckets) return -1;
|
||
for(int i=0;i<ix->n;i++){
|
||
RNbhd* r=&ix->nb[i];
|
||
for(int m=0;m<r->n_members;m++){
|
||
size_t b=geo_djb2(r->member_ids[m])%ix->nbuckets;
|
||
RE* e=malloc(sizeof*e); if(!e) continue;
|
||
e->id=r->member_ids[m]; e->nbhd=i; e->w=r->member_w[m]; e->next=ix->buckets[b]; ix->buckets[b]=e;
|
||
}
|
||
/* centered, normalized centroid for the nearest-fallback */
|
||
if(r->centroid_raw && ix->mean && ix->mean_dim==r->dim){
|
||
r->centroid_unit=malloc((size_t)r->dim*sizeof(float));
|
||
if(r->centroid_unit){
|
||
double nrm=0; for(int d=0;d<r->dim;d++){ double v=(double)r->centroid_raw[d]-ix->mean[d]; r->centroid_unit[d]=(float)v; nrm+=v*v; }
|
||
nrm=sqrt(nrm);
|
||
if(nrm>1e-12){ for(int d=0;d<r->dim;d++) r->centroid_unit[d]=(float)(r->centroid_unit[d]/nrm); }
|
||
else { free(r->centroid_unit); r->centroid_unit=NULL; }
|
||
}
|
||
}
|
||
/* fill the borrowed view */
|
||
r->view.id=r->id; r->view.hub_id=r->hub_id; r->view.n_members=r->n_members;
|
||
r->view.member_ids=r->member_ids; r->view.member_w=r->member_w;
|
||
r->view.radius=r->radius; r->view.co_registration=r->co_reg;
|
||
r->view.k_core=r->k_core; r->view.n_embedded=r->n_embedded;
|
||
}
|
||
ix->score=ix->n?calloc((size_t)ix->n,sizeof(double)):NULL;
|
||
return 0;
|
||
}
|
||
|
||
static void geo__reify_load_cb(const StoreNode* n, void* ctx){
|
||
engram_geo_reify_index_add((GeoReifyIndex*)ctx, n);
|
||
}
|
||
GeoReifyIndex* engram_geo_reify_load(EngramPagedStore* store){
|
||
if(!store) return NULL;
|
||
GeoReifyIndex* ix=engram_geo_reify_index_new(); if(!ix) return NULL;
|
||
store_scan_nodes(store, geo__reify_load_cb, ix);
|
||
if(ix->n==0 && ix->mean==NULL){ engram_geo_reify_index_free(ix); return NULL; }
|
||
engram_geo_reify_index_finalize(ix);
|
||
return ix;
|
||
}
|
||
|
||
const GeoNeighborhood* engram_geo_reify_lookup(
|
||
const GeoReifyIndex* ix,
|
||
const char* const* seed_ids, size_t n_seeds,
|
||
const float* q_emb, int q_dim){
|
||
if(!ix||ix->n<=0) return NULL;
|
||
/* (a) membership route: score each neighborhood by summed seed membership. */
|
||
if(ix->score && ix->buckets && seed_ids && n_seeds>0){
|
||
for(int i=0;i<ix->n;i++) ((GeoReifyIndex*)ix)->score[i]=0.0;
|
||
int any=0;
|
||
for(size_t s=0;s<n_seeds;s++){
|
||
const char* id=seed_ids[s]; if(!id) continue;
|
||
for(RE* e=ix->buckets[geo_djb2(id)%ix->nbuckets]; e; e=e->next)
|
||
if(strcmp(e->id,id)==0){ ((GeoReifyIndex*)ix)->score[e->nbhd]+=e->w; any=1; }
|
||
}
|
||
if(any){
|
||
int best=-1; double bv=-1;
|
||
for(int i=0;i<ix->n;i++) if(ix->score[i]>bv){ bv=ix->score[i]; best=i; }
|
||
if(best>=0 && bv>0) return &ix->nb[best].view;
|
||
}
|
||
}
|
||
/* (b) centroid-nearest fallback (centered query vs centered centroids). */
|
||
if(q_emb && q_dim>0 && ix->mean && ix->mean_dim==q_dim){
|
||
double nq=0; float* cq=malloc((size_t)q_dim*sizeof(float));
|
||
if(!cq) return NULL;
|
||
for(int d=0;d<q_dim;d++){ double v=(double)q_emb[d]-ix->mean[d]; cq[d]=(float)v; nq+=v*v; }
|
||
nq=sqrt(nq);
|
||
if(nq>1e-12){
|
||
int best=-1; double bc=-1e9;
|
||
for(int i=0;i<ix->n;i++){ RNbhd* r=&ix->nb[i]; if(!r->centroid_unit) continue;
|
||
double s=0; for(int d=0;d<q_dim;d++) s+=(double)cq[d]*r->centroid_unit[d];
|
||
s/=nq; if(s>bc){ bc=s; best=i; } }
|
||
free(cq);
|
||
if(best>=0) return &ix->nb[best].view;
|
||
} else free(cq);
|
||
}
|
||
return NULL;
|
||
}
|
||
|
||
int engram_geo_reify_count(const GeoReifyIndex* ix){ return ix?ix->n:0; }
|
||
const float* engram_geo_reify_mean(const GeoReifyIndex* ix, int* dim){
|
||
if(!ix||!ix->mean){ if(dim)*dim=0; return NULL; }
|
||
if(dim)*dim=ix->mean_dim; return ix->mean;
|
||
}
|
||
|
||
void engram_geo_reify_index_free(GeoReifyIndex* ix){
|
||
if(!ix) return;
|
||
if(ix->buckets){
|
||
for(size_t b=0;b<ix->nbuckets;b++){ RE* e=ix->buckets[b]; while(e){ RE* x=e->next; free(e); e=x; } }
|
||
free(ix->buckets);
|
||
}
|
||
for(int i=0;i<ix->n;i++){ RNbhd* r=&ix->nb[i];
|
||
free(r->id); free(r->hub_id);
|
||
for(int m=0;m<r->n_members;m++) free(r->member_ids[m]);
|
||
free(r->member_ids); free(r->member_w);
|
||
free(r->centroid_raw); free(r->centroid_unit);
|
||
}
|
||
free(ix->nb); free(ix->score); free(ix->mean);
|
||
free(ix);
|
||
}
|