/* mom.c — method-of-moments recovery: the DIRECTED read that crosses the valley.
*
* Nine local/online/greedy/SOC/wake-sleep mechanisms could not learn Mess3: the
* correct model sits across a fitness valley (partial breaks predict worse than
* uniform) and blind local search never crosses it. The resolution (evolution is
* DIRECTED, not random; learning smooths the valley -- Hinton-Nowlan): use a
* DIRECTED signal that reads structure straight from the data's global statistics,
* with NO search.
*
* Spectral / observable-operator method-of-moments (Jaeger OOM / Hsu-Kakade-Zhang):
* moments P1[i]=P(x1=i), P21[i][j]=P(x2=i,x1=j), P3x[x][i][j]=P(x3=i,x2=x,x1=j)
* operators b1 = P1 ; binf = (P21^T)^-1 P1 ; B[x] = P3x[x] * P21^-1
* predict P(a|past) ~ binf^T B[a] state ; state <- B[x] state / (binf^T B[x] state)
* With 3 states = 3 symbols and P21 full-rank, this is one 3x3 inverse. One shot.
*
* Success = log-loss -> optimal 0.794 AND the belief-geometry fractal recovered,
* in a single non-iterative directed read -- crossing what search could not.
*
* build: cc -O2 -Wall -std=c11 -o mom mom.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define K 3
static const double TT[K][K][K] = {
{{0.765,0.00375,0.00375},{0.0425,0.0675,0.00375},{0.0425,0.00375,0.0675}},
{{0.0675,0.0425,0.00375},{0.00375,0.765,0.00375},{0.00375,0.0425,0.0675}},
{{0.0675,0.00375,0.0425},{0.00375,0.0675,0.0425},{0.00375,0.00375,0.765}},
};
static unsigned long long RNG;
static unsigned long long xs(void){ RNG^=RNG<<13; RNG^=RNG>>7; RNG^=RNG<<17; return RNG; }
static double urand(void){ return (double)(xs()>>11)*(1.0/9007199254740992.0); }
static void gen_mess3(unsigned char *o,int n){ int s=0;
for(int i=0;i<n;i++){ double u=urand(),a=0; int es=0,nj=0,d=0;
for(int y=0;y<K&&!d;y++)for(int j=0;j<K&&!d;j++){a+=TT[y][s][j];if(u<a){es=y;nj=j;d=1;}}
o[i]=(unsigned char)es; s=nj; } }
static int inv3(const double A[3][3],double R[3][3]){
double det=A[0][0]*(A[1][1]*A[2][2]-A[1][2]*A[2][1])
-A[0][1]*(A[1][0]*A[2][2]-A[1][2]*A[2][0])
+A[0][2]*(A[1][0]*A[2][1]-A[1][1]*A[2][0]);
if(fabs(det)<1e-14) return 0; double id=1.0/det;
R[0][0]= (A[1][1]*A[2][2]-A[1][2]*A[2][1])*id; R[0][1]=-(A[0][1]*A[2][2]-A[0][2]*A[2][1])*id; R[0][2]= (A[0][1]*A[1][2]-A[0][2]*A[1][1])*id;
R[1][0]=-(A[1][0]*A[2][2]-A[1][2]*A[2][0])*id; R[1][1]= (A[0][0]*A[2][2]-A[0][2]*A[2][0])*id; R[1][2]=-(A[0][0]*A[1][2]-A[0][2]*A[1][0])*id;
R[2][0]= (A[1][0]*A[2][1]-A[1][1]*A[2][0])*id; R[2][1]=-(A[0][0]*A[2][1]-A[0][1]*A[2][0])*id; R[2][2]= (A[0][0]*A[1][1]-A[0][1]*A[1][0])*id;
return 1;
}
/* TRUE-model predictive distribution filter: log-loss + predictive-dist fractal */
static void true_eval(const unsigned char *s,int n,double *ll_out,long dcount[5],const int *Gs){
double b[K]={1.0/3,1.0/3,1.0/3},ll=0; long k=0;
char *seen[5]; for(int g=0;g<5;g++) seen[g]=calloc((size_t)(Gs[g]+1)*(Gs[g]+1),1);
for(int g=0;g<5;g++) dcount[g]=0;
for(int t=0;t<n;t++){ double p[K],px=0,nb[K]={0,0,0};
for(int a=0;a<K;a++){ p[a]=0; for(int i=0;i<K;i++)for(int j=0;j<K;j++){ double m=b[i]*TT[a][i][j]; p[a]+=m; nb[j]+=(a==s[t])?m:0; } px+=p[a]; }
/* predictive-dist fractal on (p0,p1) */
for(int g=0;g<5;g++){ int G=Gs[g],q0=(int)(p[0]/px*G+0.5),q1=(int)(p[1]/px*G+0.5);
if(q0<0){q0=0;} if(q0>G){q0=G;} if(q1<0){q1=0;} if(q1>G){q1=G;}
size_t key=(size_t)q0*(G+1)+q1; if(!seen[g][key]){seen[g][key]=1;dcount[g]++;} }
int x=s[t]; double pxo=p[x]>1e-300?p[x]:1e-300; ll+=-log(pxo); k++;
for(int j=0;j<K;j++) b[j]=nb[j]/pxo; /* condition on the OBSERVED symbol, not the sum */
}
for(int g=0;g<5;g++) free(seen[g]);
*ll_out=ll/(double)k;
}
int main(void){
const int NTR=6000000, NTE=2000000;
unsigned char *tr=malloc(NTR),*te=malloc(NTE);
RNG=0x9e3779b97f4a7c15ULL; gen_mess3(tr,NTR);
RNG=0xd1b54a32d192ed03ULL; gen_mess3(te,NTE);
/* ---- estimate moments from the stream (the global co-occurrence statistics) ---- */
double P1[K]={0}, P21[K][K]={{0}}, P3[K][K][K]={{{0}}};
for(int t=0;t<NTR;t++) P1[tr[t]]+=1.0;
for(int t=0;t+1<NTR;t++) P21[tr[t+1]][tr[t]]+=1.0; /* [x2][x1] */
for(int t=0;t+2<NTR;t++) P3[tr[t+1]][tr[t+2]][tr[t]]+=1.0; /* [x2][x3][x1] */
double z1=0; for(int i=0;i<K;i++) z1+=P1[i]; for(int i=0;i<K;i++) P1[i]/=z1;
double z2=0; for(int i=0;i<K;i++)for(int j=0;j<K;j++) z2+=P21[i][j];
for(int i=0;i<K;i++)for(int j=0;j<K;j++) P21[i][j]/=z2;
double z3=0; for(int x=0;x<K;x++)for(int i=0;i<K;i++)for(int j=0;j<K;j++) z3+=P3[x][i][j];
for(int x=0;x<K;x++)for(int i=0;i<K;i++)for(int j=0;j<K;j++) P3[x][i][j]/=z3;
/* ---- directed algebraic read: observable operators (one 3x3 inverse) ---- */
double P21i[3][3], P21Ti[3][3];
if(!inv3(P21,P21i)){ printf("P21 singular — not observable\n"); return 1; }
double P21T[3][3]; for(int i=0;i<3;i++)for(int j=0;j<3;j++) P21T[i][j]=P21[j][i];
inv3(P21T,P21Ti);
double b1[K]; for(int i=0;i<K;i++) b1[i]=P1[i];
double binf[K]; for(int i=0;i<K;i++){ binf[i]=0; for(int j=0;j<K;j++) binf[i]+=P21Ti[i][j]*P1[j]; }
double B[K][3][3];
for(int x=0;x<K;x++) for(int i=0;i<3;i++) for(int j=0;j<3;j++){ B[x][i][j]=0; for(int k=0;k<3;k++) B[x][i][j]+=P3[x][i][k]*P21i[k][j]; }
/* ---- evaluate: log-loss + predictive-dist fractal on the operator recursion ---- */
int Gs[5]={30,100,300,1000,3000};
char *seen[5]; for(int g=0;g<5;g++) seen[g]=calloc((size_t)(Gs[g]+1)*(Gs[g]+1),1);
long dmom[5]={0}; double st[K]; for(int i=0;i<K;i++) st[i]=b1[i];
double ll=0; long k=0;
for(int t=0;t<NTE;t++){
double num[K],tot=0;
for(int a=0;a<K;a++){ double v=0; for(int i=0;i<K;i++){ double Bi=0; for(int j=0;j<K;j++) Bi+=B[a][i][j]*st[j]; v+=binf[i]*Bi; } num[a]=v>1e-12?v:1e-12; tot+=num[a]; }
for(int g=0;g<5;g++){ int G=Gs[g],q0=(int)(num[0]/tot*G+0.5),q1=(int)(num[1]/tot*G+0.5);
if(q0<0){q0=0;} if(q0>G){q0=G;} if(q1<0){q1=0;} if(q1>G){q1=G;}
size_t key=(size_t)q0*(G+1)+q1; if(!seen[g][key]){seen[g][key]=1;dmom[g]++;} }
int x=te[t]; ll+=-log(num[x]/tot); k++;
/* state update: st <- B[x] st / (binf^T B[x] st) */
double nst[K],nrm=0; for(int i=0;i<K;i++){ nst[i]=0; for(int j=0;j<K;j++) nst[i]+=B[x][i][j]*st[j]; }
for(int i=0;i<K;i++) nrm+=binf[i]*nst[i]; if(fabs(nrm)<1e-12) nrm=1e-12;
for(int i=0;i<K;i++) st[i]=nst[i]/nrm;
}
double ll_mom=ll/(double)k;
for(int g=0;g<5;g++) free(seen[g]);
double ll_true; long dtrue[5]; true_eval(te,NTE,&ll_true,dtrue,Gs);
printf("=== METHOD-OF-MOMENTS (directed spectral read — ONE shot, no search) ===\n");
printf(" log-loss MoM (directed read) : %.4f nats\n",ll_mom);
printf(" log-loss TRUE (optimal) : %.4f nats\n",ll_true);
printf(" gap MoM - optimal : %+.4f nats (nine SEARCH methods: +0.20 .. +1.87, never 0)\n",ll_mom-ll_true);
printf(" predictive-dist FRACTAL (distinct points vs resolution):\n");
printf(" %-7s %12s %12s\n","G","MoM","TRUE");
for(int g=0;g<5;g++) printf(" G=%-5d %12ld %12ld\n",Gs[g],dmom[g],dtrue[g]);
free(tr); free(te); return 0;
}