/* learn.c — can the mating tensor be LEARNED by a local prediction-error rule?
*
* Frontier test (A): given the correct number of hidden states (3), start from a
* RANDOM mating tensor and learn it online from the Mess3 symbol stream by the
* local Hebbian rule = stochastic Baum-Welch / online EM = active-inference
* B-matrix learning:
* belief b (forward filter, current T); observe x;
* transition posterior g_ij = b_i * T[x][i][j] / P(x) (Hebbian outer product)
* accumulate counts S[x][i][j] (decaying), renormalize -> new T[.][i][.]
* belief update b'_j = sum_i g_ij
* No backward pass (filtering approximation). No knowledge of the true tensor.
*
* Success = the LEARNED filter, on fresh test data, (1) approaches the optimal
* log-loss (0.7936 nats, the true-tensor filter) and (2) still reconstructs the
* fractal belief geometry (distinct-belief count keeps growing with resolution).
* That would show the given-structure half of "learn the tensor" is solid, and
* isolate STRUCTURE emergence (how many states) as the remaining open frontier.
*
* build: cc -O2 -Wall -std=c11 -o learn learn.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define K 3
/* TRUE Mess3 tensor (a=0.9,x=0.075) — used ONLY to generate data + as oracle. */
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.0; int es=0,nj=0,done=0;
for (int sym=0; sym<K && !done; sym++) for (int j=0;j<K && !done;j++){
acc += TT[sym][s][j]; if(u<acc){ es=sym; nj=j; done=1; }
}
out[i]=(unsigned char)es; s=nj;
}
}
/* filter log-loss under a given tensor T */
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);
/* --- learner state: counts S, tensor T (init random, symmetry-broken) --- */
static double S[K][K][K], T[K][K][K];
RNG=0x243f6a8885a308d3ULL;
for(int i=0;i<K;i++){
double rs=0;
for(int s=0;s<K;s++)for(int j=0;j<K;j++){ S[s][i][j]=0.1+urand(); rs+=S[s][i][j]; }
for(int s=0;s<K;s++)for(int j=0;j<K;j++) T[s][i][j]=S[s][i][j]/rs;
}
double ll_init = filter_ll((const double(*)[K][K])T, te, NTE);
/* --- online EM / Hebbian learning over the training stream --- */
double b[K]={1.0/3,1.0/3,1.0/3};
double eta0=0.05, tau=50000.0;
for(int t=0;t<NTR;t++){
int x=tr[t];
double px=0, g[K][K];
for(int i=0;i<K;i++)for(int j=0;j<K;j++){ double m=b[i]*T[x][i][j]; g[i][j]=m; px+=m; }
if(px<1e-300)px=1e-300;
for(int i=0;i<K;i++)for(int j=0;j<K;j++) g[i][j]/=px; /* transition posterior */
double eta = eta0 * tau/(tau+t);
/* decay all counts, add the observed-symbol slice (Hebbian) */
for(int s=0;s<K;s++)for(int i=0;i<K;i++)for(int j=0;j<K;j++){
S[s][i][j] *= (1.0-eta);
if(s==x) S[s][i][j] += eta*g[i][j];
}
/* re-estimate T per source row i, with a Dirichlet/uniform floor
* (active-inference concentration prior) to prevent state collapse. */
double lam = LAM;
for(int i=0;i<K;i++){
double rs=0; for(int s=0;s<K;s++)for(int j=0;j<K;j++) rs+=S[s][i][j];
if(rs<1e-300)rs=1e-300;
for(int s=0;s<K;s++)for(int j=0;j<K;j++)
T[s][i][j]=(1.0-lam)*(S[s][i][j]/rs) + lam*(1.0/(K*K));
}
/* belief update */
double nb[K]={0,0,0}; for(int j=0;j<K;j++){ for(int i=0;i<K;i++) nb[j]+=g[i][j]; }
for(int j=0;j<K;j++) b[j]=nb[j];
}
double ll_learn = filter_ll((const double(*)[K][K])T, te, NTE);
double ll_opt = filter_ll(TT, te, NTE);
printf("=== LEARN the Mess3 mating tensor by local prediction-error (online EM / Hebbian) ===\n");
printf(" log-loss random init tensor : %.4f nats\n", ll_init);
printf(" log-loss LEARNED tensor : %.4f nats\n", ll_learn);
printf(" log-loss TRUE (optimal) tensor: %.4f nats\n", ll_opt);
printf(" gap learned - optimal : %+.4f nats\n", ll_learn-ll_opt);
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(" learned tensor (rows=source state; may be permuted vs true):\n");
for(int s=0;s<K;s++){ printf(" Tsym%d = ",s);
for(int i=0;i<K;i++){ printf("["); for(int j=0;j<K;j++) printf("%.4f%s",T[s][i][j], j<K-1?" ":""); printf("] "); }
printf("\n"); }
free(tr); free(te);
return 0;
}