/* dscale.c — calibration, not an experiment. Measures the two numbers stress.c guessed:
* (1) |delta| from a single-entry STEP perturbation near the optimum -> what T must sit below
* (2) window sampling noise in the log-loss estimate -> the floor on resolvable delta
* If (2) >= (1), the acceptance test is a coin flip and no temperature schedule can work. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.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(unsigned char*o,int n){int s=0;for(int i=0;i<n;i++){double u=urand(),a=0;int e=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){e=y;nj=j;d=1;}}o[i]=(unsigned char)e;s=nj;}}
static double ll_(const double T[K][K][K],const unsigned char*s,long off,int n){
double b[K]={1.0/3,1.0/3,1.0/3},L=0;
for(int t=0;t<n;t++){int x=s[off+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;
L+=-log(px);
for(int j=0;j<K;j++)b[j]=nb[j]/px;}
return L/(double)n;}
static void rn(double T[K][K][K]){for(int i=0;i<K;i++){double z=0;
for(int s=0;s<K;s++)for(int j=0;j<K;j++){if(T[s][i][j]<1e-6)T[s][i][j]=1e-6;z+=T[s][i][j];}
for(int s=0;s<K;s++)for(int j=0;j<K;j++)T[s][i][j]/=z;}}
int main(void){
const long NTR=600000;unsigned char*tr=malloc(NTR);
RNG=0x9e3779b97f4a7c15ULL;gen(tr,(int)NTR);
int WINS[3]={10000,50000,200000};
RNG=0xabcdef1234567ULL;
for(int w=0;w<3;w++){
int W=WINS[w];
/* (1) |delta| from STEP=0.06 single-entry perturbation, SAME window (common random numbers) */
double s1=0,s2=0;int N=300;
for(int k=0;k<N;k++){
long off=(long)(urand()*(double)(NTR-W-1));
double base=ll_(TT,tr,off,W),T2[K][K][K];memcpy(T2,TT,sizeof T2);
int si=(int)(urand()*K),ii=(int)(urand()*K),jj=(int)(urand()*K);
if(si>=K)si=K-1;
if(ii>=K)ii=K-1;
if(jj>=K)jj=K-1;
T2[si][ii][jj]*=exp(0.06*(2*urand()-1));rn(T2);
double d=fabs(ll_((const double(*)[K][K])T2,tr,off,W)-base);
s1+=d;s2+=d*d;}
double md=s1/N,sd=sqrt(s2/N-md*md);
/* (2) window-to-window noise in the SAME tensor's score */
double n1=0,n2=0;int NN=300;
for(int k=0;k<NN;k++){long off=(long)(urand()*(double)(NTR-W-1));
double v=ll_(TT,tr,off,W);n1+=v;n2+=v*v;}
double mn=n1/NN,sn=sqrt(n2/NN-mn*mn);
printf("WIN=%-7d |delta| from one STEP: %.5f (sd %.5f) window noise sd: %.5f ratio noise/delta = %.2f\n",
W,md,sd,sn,sn/md);
}
free(tr);return 0;}