/* sns.c — the nonunifilar falsifier.
*
* Runs our finite prediction-suffix-tree method (the same split+merge as
* stream.c) against the SIMPLE NONUNIFILAR SOURCE (SNS), the canonical process
* with INFINITELY many causal states, and compares it head-to-head with the
* EXACT belief-state filter (the mixed-state presentation / optimal predictor).
*
* SNS (hidden states A,B; start A):
* A: emit 0 -> A (p=1/2) | emit 1 -> B (p=1/2)
* B: emit 1 -> A (p=1/2) | emit 1 -> B (p=1/2) (B always emits 1)
* Nonunifilar: from B, symbol 1 leads to A OR B. A '0' pins you to A (reset);
* so the causal state = length of the current run of 1s since the last 0 ->
* a COUNTABLY INFINITE set of states, predictions accumulating to a fixed point.
*
* Prediction (the falsification target): a finite-depth PST CANNOT hold an
* infinite causal-state set. It must truncate at DMAX and lump every longer run
* together. The exact belief filter has no such ceiling — it is a 1-D dynamical
* system on the belief simplex whose visited points are the (infinite) causal
* states.
*
* build: cc -O2 -Wall -std=c11 -o sns sns.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DMAX 20
#define MINOBS 50
#define CHI2_CRIT 10.828
#define ALPHA 1.0
/* ---- deterministic PRNG ---- */
static unsigned long long RNG;
static unsigned long long xs(void){ RNG^=RNG<<13; RNG^=RNG>>7; RNG^=RNG<<17; return RNG; }
static int coin(void){ return (int)(xs() & 1); }
/* ---- SNS generator ---- */
static void gen_sns(unsigned char *b, int n){
int s = 0; /* 0 = A, 1 = B; start in A */
for (int i = 0; i < n; i++){
if (s == 0){ /* A */
if (coin()){ b[i] = 1; s = 1; } /* emit 1 -> B */
else { b[i] = 0; s = 0; } /* emit 0 -> A */
} else { /* B: always emit 1, jump A/B */
b[i] = 1; s = coin();
}
}
}
/* ============ finite PST method (identical to stream.c) ============ */
typedef struct Node { long c[2]; struct Node *ch[2]; int minted; } Node;
static Node *newnode(void){ return calloc(1, sizeof(Node)); }
static long total(const Node *n){ return n->c[0] + n->c[1]; }
static double p1(const Node *n){ return (n->c[1]+ALPHA) / (total(n) + 2*ALPHA); }
static double pof(const Node *n, int s){ return s ? p1(n) : 1.0 - p1(n); }
static void build(Node *root, const unsigned char *s, int n){
for (int t = 1; t < n; t++){
Node *node = root; node->c[s[t]]++;
for (int L = 1; L <= DMAX && t - L >= 0; L++){
int a = s[t-L];
if (!node->ch[a]) node->ch[a] = newnode();
node = node->ch[a];
node->c[s[t]]++;
}
}
}
static double gtest(const Node *child, const Node *parent){
long N = total(child); if (!N) return 0.0;
double G = 0.0;
for (int s = 0; s < 2; s++)
if (child->c[s] > 0) G += (double)child->c[s] * log((double)child->c[s] / ((double)N * pof(parent, s)));
return 2.0 * G;
}
static int grow(Node *root){
root->minted = 1;
Node *queue[1<<16]; int head=0, tail=0; queue[tail++]=root;
int created = 0;
while (head < tail){
Node *node = queue[head++];
for (int a = 0; a < 2; a++){
Node *ch = node->ch[a];
if (ch && total(ch) >= MINOBS && gtest(ch, node) >= CHI2_CRIT){
ch->minted = 1; created++;
if (tail < (1<<16)) queue[tail++] = ch;
}
}
}
return created;
}
static const Node *ctx(const Node *root, const unsigned char *s, int t){
const Node *node = root;
for (int L = 1; L <= DMAX && t - L >= 0; L++){
int a = s[t-L];
if (node->ch[a] && node->ch[a]->minted) node = node->ch[a];
else break;
}
return node;
}
static double logloss_pst(const Node *root, const unsigned char *s, int n){
double ll = 0; long k = 0;
for (int t = DMAX; t < n; t++){ ll += -log(pof(ctx(root, s, t), s[t])); k++; }
return ll / (double)k;
}
static long LC_[4096][2]; static int NL_;
static void collect(Node *n){
if(!n) return;
if(n->minted){
int deeper=(n->ch[0]&&n->ch[0]->minted)||(n->ch[1]&&n->ch[1]->minted);
if(!deeper && NL_<4096){ LC_[NL_][0]=n->c[0]; LC_[NL_][1]=n->c[1]; NL_++; }
}
collect(n->ch[0]); collect(n->ch[1]);
}
static double homog_G(long a0,long a1,long b0,long b1){
double N=(double)(a0+a1+b0+b1); if(N==0) return 0;
double rA=a0+a1, rB=b0+b1, c0=a0+b0, c1=a1+b1, G=0;
if(a0>0){ double E=rA*c0/N; G += a0*log(a0/E); }
if(a1>0){ double E=rA*c1/N; G += a1*log(a1/E); }
if(b0>0){ double E=rB*c0/N; G += b0*log(b0/E); }
if(b1>0){ double E=rB*c1/N; G += b1*log(b1/E); }
return 2*G;
}
static int minimize(Node *root){
NL_=0; collect(root);
int alive[4096]; for(int i=0;i<NL_;i++) alive[i]=1;
int nstates=NL_;
for(;;){
double best=1e18; int bi=-1,bj=-1;
for(int i=0;i<NL_;i++) if(alive[i])
for(int j=i+1;j<NL_;j++) if(alive[j]){
double g=homog_G(LC_[i][0],LC_[i][1],LC_[j][0],LC_[j][1]);
if(g<best){ best=g; bi=i; bj=j; }
}
if(bi<0 || best>=CHI2_CRIT) break;
LC_[bi][0]+=LC_[bj][0]; LC_[bi][1]+=LC_[bj][1]; alive[bj]=0; nstates--;
}
return nstates;
}
/* ============ EXACT belief-state filter (the MSP / optimal predictor) ============
* belief bA = P(state==A | history). A '0' resets bA->1. A '1' updates by Bayes.
* P(1|b) = bA/2 + bB ; P(0|b) = bA/2.
* after 1: bA' = (bB/2)/P(1) ; after 0: bA' = 1.
*/
static double filter_logloss(const unsigned char *s, int n, double *out_states, long *out_distinct){
double bA = 1.0; /* start in A */
double ll = 0; long k = 0;
/* quantize visited beliefs to count distinct causal states (1e-4 grid). */
static int seen[100001];
for (int i=0;i<=100000;i++) seen[i]=0;
long distinct = 0;
for (int t = 0; t < n; t++){
double bB = 1.0 - bA;
double p1b = bA*0.5 + bB; /* P(next=1) */
int q = (int)(bA*100000.0 + 0.5); if(q<0)q=0; if(q>100000)q=100000;
if(!seen[q]){ seen[q]=1; distinct++; }
int x = s[t];
double px = x ? p1b : (bA*0.5);
if (px < 1e-15) px = 1e-15;
ll += -log(px); k++;
/* update */
if (x == 0) bA = 1.0;
else bA = (bB*0.5) / p1b;
}
*out_distinct = distinct;
(void)out_states;
return ll / (double)k;
}
int main(void){
const int NTR = 2000000, NTE = 400000;
unsigned char *tr = malloc(NTR), *te = malloc(NTE);
RNG = 0x9e3779b97f4a7c15ULL; gen_sns(tr, NTR);
RNG = 0xd1b54a32d192ed03ULL; gen_sns(te, NTE);
/* ---- our finite PST method ---- */
Node *root = newnode();
build(root, tr, NTR);
int created = grow(root);
int nmin = minimize(root);
double ll0 = 0; long k=0;
for (int t=0;t<NTE;t++){ ll0 += -log(pof(root, te[t])); k++; } ll0/=(double)k;
double llp = logloss_pst(root, te, NTE);
/* ---- exact belief filter ---- */
long distinct=0;
double llf = filter_logloss(te, NTE, NULL, &distinct);
printf("=== SIMPLE NONUNIFILAR SOURCE (infinite causal states) ===\n");
printf(" finite PST (our stream.c method):\n");
printf(" distinctions created (split) : %d\n", created);
printf(" minimized to causal states : %d (DMAX=%d ceiling)\n", nmin, DMAX);
printf(" log-loss order-0 baseline : %.4f nats\n", ll0);
printf(" log-loss grown world-model : %.4f nats\n", llp);
printf(" exact belief filter (MSP / optimal predictor):\n");
printf(" log-loss (= entropy rate est) : %.4f nats\n", llf);
printf(" distinct belief states visited : %ld (the infinite/accumulating causal-state set)\n", distinct);
printf(" gap PST - optimal : %+.4f nats\n", llp - llf);
/* accumulation table: P(next=1) after a run of k ones since the last 0. */
printf(" belief accumulation P(next=1 | run of k ones since last 0):\n");
double bA = 1.0; /* just saw a 0 -> state A */
for (int kk = 0; kk <= 8; kk++){
double bB = 1.0 - bA;
double p = bA*0.5 + bB;
printf(" k=%d : P(1)=%.5f bA=%.5f\n", kk, p, bA);
bA = (bB*0.5) / p; /* advance one more '1' */
}
printf(" ... (never exactly repeats -> countably INFINITE causal states)\n\n");
free(tr); free(te);
return 0;
}