/* wakesleep2.c — the consolidation ratchet on top of the SOC controller.
*
* SOC proved: at beta_c the system VISITS the correct broken state but online
* averaging washes it out (fluctuations are symmetric, they average away). The
* ratchet adds hysteresis so a correct break is FROZEN instead of averaged:
*
* WAKE : run online DA-EM at SOC-controlled beta (explores, breaks symmetry
* transiently — the "mutation" source = critical fluctuations).
* SLEEP : every EVAL_EVERY steps, evaluate the current model's log-loss on a
* held-out batch. If it IMPROVED on the best-so-far, CONSOLIDATE
* (best = current). If it got WORSE, ROLL BACK (current = best). The
* rollback is the anti-washout: the model can ratchet UP when a critical
* fluctuation breaks correctly, but can never wash back DOWN.
*
* This is wake(mutate)/sleep(select) — an evolutionary ratchet driven by SOC
* critical fluctuations. Success = best-model log-loss -> optimal 0.794 AND the
* true fractal, which pure online (soc.c) could not hold.
*
* build: cc -O2 -Wall -std=c11 -o wakesleep2 wakesleep.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define K 3
#ifndef ETA
#define ETA 0.001
#endif
#ifndef KAPPA
#define KAPPA 0.0005
#endif
#ifndef NOISE
#define NOISE 0.02
#endif
#ifndef WIN
#define WIN 50000
#endif
#ifndef STEP
#define STEP 0.02
#endif
#ifndef EVAL_EVERY
#define EVAL_EVERY 20000
#endif
#define EVAL_N 100000
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; }
/* SLEEP = batch-EM replay to convergence: fully DESCEND whatever partial break
* wake seeded, to its nearest fixed point, before judging it (crosses the valley). */
static void batch_em(double T[K][K][K], const unsigned char *buf, int n, int iters){
static double C[K][K][K];
for(int it=0;it<iters;it++){
memset(C,0,sizeof C);
double b[K]={1.0/3,1.0/3,1.0/3};
for(int t=0;t<n;t++){ int x=buf[t]; 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;}
for(int i=0;i<K;i++)for(int j=0;j<K;j++) C[x][i][j]+=m[i][j]/px;
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; } }
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+=C[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]=(C[s][i][j]+KAPPA)/den; } }
}
int main(void){
const int NTR=6000000;
unsigned char *tr=malloc(NTR),*ev=malloc(EVAL_N),*te=malloc(2000000);
RNG=0x9e3779b97f4a7c15ULL; gen_mess3(tr,NTR);
RNG=0x5555555555555555ULL; gen_mess3(ev,EVAL_N); /* held-out sleep-eval batch */
RNG=0xd1b54a32d192ed03ULL; gen_mess3(te,2000000); /* final test */
static double N[K][K][K],T[K][K][K],Nbest[K][K][K],Tbest[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];} }
memcpy(Nbest,N,sizeof N); memcpy(Tbest,T,sizeof T);
double Lbest=filter_ll((const double(*)[K][K])T,ev,EVAL_N);
double Linit=Lbest;
double b[K]={1.0/3,1.0/3,1.0/3};
double beta=0.90,dir=1.0,chi_prev=-1.0,ws=0,ws2=0; long wc=0;
int accepts=0,rollbacks=0;
for(int t=0;t<NTR;t++){
int x=tr[t];
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;}
double r[K]; for(int i=0;i<K;i++){ double si=0; for(int j=0;j<K;j++) si+=m[i][j]; r[i]=si/px; }
double s_op=0; for(int i=0;i<K;i++){ double d=r[i]-1.0/K; s_op+=d*d; }
ws+=s_op; ws2+=s_op*s_op; wc++;
if(wc>=WIN){ double mean=ws/wc,chi=ws2/wc-mean*mean;
if(chi_prev>=0&&chi<chi_prev) dir=-dir;
beta+=dir*STEP; if(beta<0.5)beta=0.5; if(beta>2.0)beta=2.0;
chi_prev=chi; ws=0; ws2=0; wc=0; }
/* WAKE: online DA-EM at SOC beta */
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;
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; }
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; }
/* SLEEP: batch-EM replay-to-convergence (consolidate), THEN keep-best/rollback */
if(t>0 && t%EVAL_EVERY==0){
batch_em(T, ev, EVAL_N, 20); /* descend the wake perturbation */
double L=filter_ll((const double(*)[K][K])T,ev,EVAL_N);
if(L<Lbest-1e-6){ Lbest=L; memcpy(Tbest,T,sizeof T); accepts++; }
else { memcpy(T,Tbest,sizeof T); rollbacks++; }
memcpy(N,T,sizeof N); /* wake resumes from consolidated model */
if(t%600000==0){ double diag=(Tbest[0][0][0]+Tbest[1][1][1]+Tbest[2][2][2])/3.0;
fprintf(stderr," t=%7d beta=%.2f L_best=%.4f diag_best=%.3f (acc=%d roll=%d)\n",
t,beta,Lbest,diag,accepts,rollbacks); } }
}
double ll_learn=filter_ll((const double(*)[K][K])Tbest,te,2000000);
double ll_opt =filter_ll(TT,te,2000000);
printf("=== WAKE/SLEEP CONSOLIDATION RATCHET on Mess3 ===\n");
printf(" log-loss init random : %.4f nats\n",Linit);
printf(" log-loss BEST (ratcheted) : %.4f nats\n",ll_learn);
printf(" log-loss TRUE (optimal) : %.4f nats\n",ll_opt);
printf(" gap best - optimal : %+.4f nats (soc/online washout was +0.30)\n",ll_learn-ll_opt);
printf(" ratchet: %d consolidations, %d rollbacks\n",accepts,rollbacks);
double diag=(Tbest[0][0][0]+Tbest[1][1][1]+Tbest[2][2][2])/3.0;
printf(" best-model diagonal : %.3f\n",diag);
printf(" FRACTAL PROBE (best model) 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])Tbest,te,2000000,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(ev); free(te); return 0;
}