/* da-em.c — the synthesis: deterministic-annealing EM (data-fit + symmetry-breaking).
*
* The triangulation:
* EM alone (learn.c/critical.c) -> washout/collapse : data-fit, NO symmetry-breaking.
* competition alone (competitive.c) -> confident-wrong : symmetry-breaking, NO data-fit anchor.
* => need BOTH, balanced at criticality. That balance is TEMPERATURE, annealed
* slowly through the bifurcation (Ueda-Nakano deterministic-annealing EM).
*
* inference : soft Bayesian belief filter (untempered) -> traces the fractal.
* M-step : EM soft counts (the DATA-FIT anchor) with responsibilities
* g_ij ~ m_ij^beta (tempered). beta LOW = high temp = symmetric basin;
* as beta crosses beta_c the symmetric solution BIFURCATES and states
* specialize -- but anchored by the likelihood, so onto the CORRECT branch.
* schedule : beta annealed SLOWLY beta0 -> beta1 (not slammed like the homeostat),
* constant eta (learning doesn't freeze before annealing ends),
* tiny persistent NOISE = the symmetry-breaking seed to amplify.
*
* Success target (BOTH, since competitive.c faked the geometry with a wrong tensor):
* log-loss -> optimal 0.794 AND fractal recovered (distinct-belief count grows).
*
* build: cc -O2 -Wall -std=c11 -o da-em da-em.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define K 3
#ifndef BETA0
#define BETA0 0.6
#endif
#ifndef BETA1
#define BETA1 1.2
#endif
#ifndef ETA
#define ETA 0.001
#endif
#ifndef KAPPA
#define KAPPA 0.0005
#endif
#ifndef NOISE
#define NOISE 0.02
#endif
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 *out,int n){
int s=0;
for(int i=0;i<n;i++){ double u=urand(),acc=0; int es=0,nj=0,d=0;
for(int sym=0;sym<K&&!d;sym++)for(int j=0;j<K&&!d;j++){acc+=TT[sym][s][j];if(u<acc){es=sym;nj=j;d=1;}}
out[i]=(unsigned char)es; s=nj; }
}
static double filter_ll(const double T[K][K][K],const unsigned char *s,int n){
double b[K]={1.0/3,1.0/3,1.0/3},ll=0; long k=0;
for(int t=0;t<n;t++){ int x=s[t]; double px=0,nb[K]={0,0,0};
for(int i=0;i<K;i++)for(int j=0;j<K;j++){double m=b[i]*T[x][i][j];nb[j]+=m;px+=m;}
if(px<1e-300){px=1e-300;}
ll+=-log(px);k++;
for(int j=0;j<K;j++)b[j]=nb[j]/px; }
return ll/(double)k;
}
static long filter_distinct(const double T[K][K][K],const unsigned char *s,int n,int G){
double b[K]={1.0/3,1.0/3,1.0/3}; char *seen=calloc((size_t)(G+1)*(G+1),1); long d=0;
for(int t=0;t<n;t++){ int q0=(int)(b[0]*G+0.5),q1=(int)(b[1]*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[key]){seen[key]=1;d++;}
int x=s[t]; double px=0,nb[K]={0,0,0};
for(int i=0;i<K;i++)for(int j=0;j<K;j++){double m=b[i]*T[x][i][j];nb[j]+=m;px+=m;}
if(px<1e-300){px=1e-300;}
for(int j=0;j<K;j++)b[j]=nb[j]/px; }
free(seen); return d;
}
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);
static double N[K][K][K],T[K][K][K];
RNG=0x243f6a8885a308d3ULL;
for(int i=0;i<K;i++){ double z=0;
for(int s=0;s<K;s++)for(int j=0;j<K;j++){N[s][i][j]=0.1+urand();z+=N[s][i][j];}
for(int s=0;s<K;s++)for(int j=0;j<K;j++){N[s][i][j]/=z; T[s][i][j]=N[s][i][j];} }
double ll_init=filter_ll((const double(*)[K][K])T,te,NTE);
double b[K]={1.0/3,1.0/3,1.0/3};
for(int t=0;t<NTR;t++){
int x=tr[t];
/* proper DA recipe (beta_c ~ 1.0, measured): OVERSHOOT to break onto the
* principal branch, then RELAX to beta=1.0 for unbiased ML refinement. */
double a=(double)t/1.0e6; /* millions of steps */
double beta;
if(a<2.0) beta = 0.7 + (1.12-0.7)*(a/2.0); /* 0.7 -> 1.12 : break */
else if(a<3.0) beta = 1.12 - (1.12-1.0)*(a-2.0); /* 1.12 -> 1.0 : relax */
else beta = 1.0; /* hold 1.0 : refine */
double m[K][K],px=0;
for(int i=0;i<K;i++)for(int j=0;j<K;j++){ m[i][j]=b[i]*T[x][i][j]; px+=m[i][j]; }
if(px<1e-300){px=1e-300;}
/* tempered EM responsibilities + tiny symmetry-breaking noise */
double g[K][K],gs=0;
for(int i=0;i<K;i++)for(int j=0;j<K;j++){ g[i][j]=pow(m[i][j],beta)*(1.0+NOISE*(urand()-0.5)); gs+=g[i][j]; }
if(gs<1e-300){gs=1e-300;}
for(int i=0;i<K;i++)for(int j=0;j<K;j++) g[i][j]/=gs;
/* decaying-average sufficient stats (constant eta, only observed symbol gains) */
for(int s=0;s<K;s++)for(int i=0;i<K;i++)for(int j=0;j<K;j++){
N[s][i][j]*=(1.0-ETA); if(s==x) N[s][i][j]+=ETA*g[i][j];
}
for(int i=0;i<K;i++){ double z=0; for(int s=0;s<K;s++)for(int j=0;j<K;j++) z+=N[s][i][j];
double den=z+KAPPA*K*K; if(den<1e-300)den=1e-300;
for(int s=0;s<K;s++)for(int j=0;j<K;j++) T[s][i][j]=(N[s][i][j]+KAPPA)/den; }
/* untempered belief update */
double nb[K]={0,0,0}; for(int j=0;j<K;j++){ for(int i=0;i<K;i++) nb[j]+=m[i][j]; b[j]=nb[j]/px; }
if(t%500000==0){
double diag=(T[0][0][0]+T[1][1][1]+T[2][2][2])/3.0;
fprintf(stderr," t=%7d beta=%.2f diag(true~.765)=%.3f b=[%.2f %.2f %.2f]\n",t,beta,diag,b[0],b[1],b[2]);
}
}
double ll_learn=filter_ll((const double(*)[K][K])T,te,NTE);
double ll_opt =filter_ll(TT,te,NTE);
printf("=== DA-EM SYNTHESIS on Mess3 (beta %.2f->%.2f, eta=%.4f, noise=%.2f) ===\n",(double)BETA0,(double)BETA1,(double)ETA,(double)NOISE);
printf(" log-loss random init : %.4f nats\n",ll_init);
printf(" log-loss LEARNED (DA-EM) : %.4f nats\n",ll_learn);
printf(" log-loss TRUE (optimal) : %.4f nats\n",ll_opt);
printf(" gap learned - optimal : %+.4f nats (collapse +0.20 ; washout +0.30 ; confident-wrong +1.87)\n",ll_learn-ll_opt);
double diag=(T[0][0][0]+T[1][1][1]+T[2][2][2])/3.0;
printf(" symmetry broken CORRECTLY? : diag %.3f (true 0.765 ; uniform 0.111 ; anti 0.000)\n",diag);
printf(" FRACTAL PROBE (learned tensor) distinct belief states vs resolution:\n");
int Gs[]={30,100,300,1000,3000}; long prev=0; int prevG=0;
for(int gi=0;gi<5;gi++){ long d=filter_distinct((const double(*)[K][K])T,te,NTE,Gs[gi]);
double dim=(prev&&prevG)?log((double)d/prev)/log((double)Gs[gi]/prevG):0.0;
printf(" G=%-5d : %8ld distinct (box-dim slope %.2f)\n",Gs[gi],d,dim); prev=d;prevG=Gs[gi]; }
printf(" [TRUE Mess3: 396/2399/8924/26394/54821]\n");
free(tr); free(te); return 0;
}