/* soc.c — the self-organized-criticality controller.
*
* DA-EM proved the correct emergent model lives only AT beta_c (~1.0), flanked by
* washout (below) and confident-wrong (above), and no FIXED or SCHEDULED beta sits
* there stably. SOC = let the learner FIND and HOLD beta_c from its own dynamics.
*
* Signal: order parameter s = sum_i (r_i - 1/K)^2 (0 = symmetric, high = broken).
* Its temporal variance chi = susceptibility, which PEAKS at the critical point.
* Controller: perturb-and-observe peak-tracker (MPPT-style) on chi -- every window,
* if chi rose keep moving beta the same way, else reverse. beta self-organizes to
* argmax_beta chi = beta_c, with NO set-point (the thing usage-entropy homeostasis
* lacked). Learning is DA-EM's tempered soft-EM, beta now driven by the controller.
*
* Honest open risk (tested, not assumed): the OVER-broken phase also fluctuates
* (belief hops between corners), so chi may not peak cleanly at beta_c. Watch whether
* beta parks near ~1.0 with log-loss -> 0.794 AND the true fractal, or drifts high.
*
* build: cc -O2 -Wall -std=c11 -o soc soc.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.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
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 b[K]={1.0/3,1.0/3,1.0/3};
/* SOC controller state */
double beta=0.90, dir=1.0, chi_prev=-1.0;
double ws=0.0, ws2=0.0; long wc=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; }
/* order parameter + susceptibility accumulation */
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; /* windowed variance = susceptibility */
if(chi_prev>=0.0 && chi<chi_prev) dir=-dir; /* susceptibility fell -> reverse */
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;
}
/* DA-EM learning at the SOC-controlled 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; }
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=%.3f chi=%.4f s_op=%.3f diag=%.3f\n",t,beta,chi_prev,s_op,diag);
}
}
double ll_learn=filter_ll((const double(*)[K][K])T,te,NTE);
double ll_opt =filter_ll(TT,te,NTE);
printf("=== SOC-VARIANCE CONTROLLER on Mess3 ===\n");
printf(" log-loss LEARNED (SOC) : %.4f nats\n",ll_learn);
printf(" log-loss TRUE (optimal) : %.4f nats\n",ll_opt);
printf(" gap learned - optimal : %+.4f nats (best prior: collapse/over-hard +0.20)\n",ll_learn-ll_opt);
printf(" final beta parked at : %.3f (measured beta_c ~ 1.0)\n",beta);
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;
}