/* stream.c — the create-a-distinction-from-error step, standalone.
*
* A basal agent watches a binary stream and grows a world model from ZERO by
* one mechanism: start with the empty context (order-0), and CREATE a new
* distinction ("the symbol at lag L was a") only when adding it *significantly*
* changes the next-symbol prediction (G-test vs the parent context). Each new
* distinction exposes finer contexts with their own residual error, which drives
* further splits — autocatalysis — halting where the residual is only noise.
*
* Test: the EVEN PROCESS (1s occur in even-length blocks). It has exactly 2
* causal states (parity of the trailing 1-run) and NO finite-order Markov model
* captures it. Control: an i.i.d. coin, which has NO structure — the agent must
* grow nothing.
*
* Grounding: this is prediction-suffix-tree / causal-state reconstruction; its
* fixed point is the minimal sufficient statistic for predicting the stream.
*
* build: cc -O2 -Wall -std=c11 -o stream stream.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DMAX 20
#define MINOBS 50
#define CHI2_CRIT 10.828 /* chi-squared df=1, p<0.001 */
#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); }
/* ---- stream generators ---- */
static void gen_even(unsigned char *b, int n){
int in_s1 = 0; /* S1 = mid-pair, must emit 1 */
for (int i = 0; i < n; i++){
int x;
if (in_s1){ x = 1; in_s1 = 0; }
else if (coin()){ x = 1; in_s1 = 1; } /* open a 1-block */
else { x = 0; } /* stay in S0 */
b[i] = (unsigned char)x;
}
}
static void gen_iid(unsigned char *b, int n){ for (int i=0;i<n;i++) b[i]=(unsigned char)coin(); }
/* ---- context trie ---- */
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); }
/* accumulate next-symbol counts for every context up to DMAX. */
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]]++;
}
}
}
/* G = 2 sum_s O_s ln(O_s / (N * P_parent(s))) — does this context predict
* differently from its parent by more than sampling noise? */
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;
}
/* the catalytic step: mint a child context only when the split is significant. */
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; /* number of distinctions created (excl. root) */
}
/* prediction: walk to the deepest MINTED context matching the history. */
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;
}
/* mean predictive log-loss (nats) on a held-out stream. */
static double logloss(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;
}
/* P(next=1 | recent history), pattern given most-recent-first. */
static double predict_after(const Node *root, const int *pat, int len){
const Node *node = root;
for (int L = 1; L <= len; L++){
int a = pat[L-1];
if (node->ch[a] && node->ch[a]->minted) node = node->ch[a];
else break;
}
return p1(node);
}
/* ---- the merge dual: collapse contexts that predict identically ---- */
/* homogeneity G-test between two count vectors (2x2 contingency, df=1).
* G < CHI2_CRIT => distributions are statistically the SAME => merge. */
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 long LC_[4096][2]; static int NL_;
static void collect(Node *n){ /* gather the minted LEAVES = prediction states */
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]);
}
/* agglomeratively merge contexts whose predictions are indistinguishable. */
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--;
}
printf(" minimized to causal states : %d [", nstates);
int first=1;
for(int i=0;i<NL_;i++) if(alive[i]){
long tot=LC_[i][0]+LC_[i][1];
printf("%sP(1)=%.2f(n=%ld)", first?"":", ", (double)LC_[i][1]/tot, tot);
first=0;
}
printf("]\n");
return nstates;
}
int main(void){
const int NTR = 500000, NTE = 100000;
unsigned char *tr = malloc(NTR), *te = malloc(NTE);
for (int which = 0; which < 2; which++){
const char *name = which ? "i.i.d. coin (control: no structure)"
: "EVEN PROCESS (2 causal states)";
RNG = 0x9e3779b97f4a7c15ULL; gen_even(tr, NTR); if (which) gen_iid(tr, NTR);
RNG = 0xd1b54a32d192ed03ULL; gen_even(te, NTE); if (which) gen_iid(te, NTE);
Node *root = newnode();
build(root, tr, NTR);
int created = grow(root);
/* order-0 baseline: predict from root only. */
double ll0 = 0; long k = 0;
for (int t = 0; t < NTE; t++){ ll0 += -log(pof(root, te[t])); k++; }
ll0 /= (double)k;
double llm = logloss(root, te, NTE);
printf("=== %s ===\n", name);
printf(" distinctions created (from zero) : %d\n", created);
minimize(root);
printf(" log-loss order-0 baseline : %.4f nats\n", ll0);
printf(" log-loss grown world-model : %.4f nats (%+.4f)\n", llm, llm - ll0);
printf(" [optimal for even process ~0.462; i.i.d ~0.693 = ln2]\n");
if (!which){
printf(" discovered parity structure P(next=1 | recent history):\n");
int a[] = {0}; printf(" after ...0 : %.3f (S0: 0/1 free)\n", predict_after(root,a,1));
int b[] = {1,0}; printf(" after ...01 : %.3f (odd run -> MUST be 1)\n", predict_after(root,b,2));
int c[] = {1,1,0}; printf(" after ...011 : %.3f (even run -> 0/1 free)\n", predict_after(root,c,3));
int d[] = {1,1,1,0}; printf(" after ...0111 : %.3f (odd run -> MUST be 1)\n", predict_after(root,d,4));
int e[] = {1,1,1,1,0}; printf(" after ...01111 : %.3f (even run -> 0/1 free)\n", predict_after(root,e,5));
}
printf("\n");
}
free(tr); free(te);
return 0;
}