M9 §5: geometry OPERATORS as C functions + EL builtins (read-only, staged)
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.
This commit is contained in:
@@ -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;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.
|
||||
* ═══════════════════════════════════════════════════════════════════════════ */
|
||||
|
||||
Reference in New Issue
Block a user