/* poscheck.c — the positive control stress.c was missing.
* Q1: can the scorer recognize the optimum (does a TRUE-seeded unit stay)?
* Q2: can this Metropolis reach it from NEAR-true (is the search capable at all)?
* If Q2 fails, stress.c's null is UNINFORMATIVE, not a null. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define K 3
#ifndef TEMP
#define TEMP 0.0050
#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(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;const int NTE=200000;
unsigned char*tr=malloc(NTR),*te=malloc(NTE);
RNG=0x9e3779b97f4a7c15ULL;gen(tr,(int)NTR); RNG=0xd1b54a32d192ed03ULL;gen(te,NTE);
double opt=ll_(TT,te,0,NTE);
printf("optimal (test) %.4f nats\n",opt);
/* Q1: TRUE-seeded, run the SAME Metropolis stress.c uses */
for(int trial=0;trial<2;trial++){
double T[K][K][K],Tp[K][K][K];
memcpy(T,TT,sizeof T);
RNG=0x243f6a8885a308d3ULL+trial;
if(trial==1){ /* NEAR-true: perturb every entry by up to +-25% */
for(int s=0;s<K;s++)for(int i=0;i<K;i++)for(int j=0;j<K;j++)T[s][i][j]*=exp(0.25*(2*urand()-1));
rn(T);}
double start=ll_((const double(*)[K][K])T,te,0,NTE);
long off0=0; double cur=ll_((const double(*)[K][K])T,tr,off0,10000);
long up=0;
for(int sw=0;sw<4000;sw++){
long off=(long)(urand()*(double)(NTR-10001));
cur=ll_((const double(*)[K][K])T,tr,off,10000);
memcpy(Tp,T,sizeof Tp);
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;
Tp[si][ii][jj]*=exp(0.06*(2*urand()-1)); rn(Tp);
double p=ll_((const double(*)[K][K])Tp,tr,off,10000),d=p-cur;
if(d<=0||urand()<exp(-d/TEMP)){if(d>0)up++;memcpy(T,Tp,sizeof T);}
}
double end=ll_((const double(*)[K][K])T,te,0,NTE);
printf("%-10s start %.4f (gap %+.4f) -> end %.4f (gap %+.4f) uphill accepts %ld/4000\n",
trial?"NEAR-true":"TRUE-seed",start,start-opt,end,end-opt,up);
}
free(tr);free(te);return 0;}