/* ribophi.c — Phi_r (causal emergence) probe on ribozyme.c's replicase bootstrap.
*
* WHY. Pigozzi & Levin (arXiv 2607.28250) report that in the GARD model, spikes in the
* PhiID causal-emergence atom Phi_r PRECEDE and predict the appearance of self-replicators,
* and that Phi_r correlates with no other network/dynamical metric they tested. Our
* ribozyme.c has a SHARP known ground truth: the replicase bootstrap IGNITES for compartment
* size 4 and 8 (phi 6 -> ~11.5) and FAILS well-mixed and at size 16 and 32 (phi stays ~5.5-5.9).
* That makes it a natural test bed: if Phi_r is a real leading indicator of an impending
* transition, it should separate those arms EARLY, before the outcome is visible.
*
* WHAT IS COPIED VERBATIM from ribozyme.c: the RNG, phi(), replicate(), and the whole run()
* dynamics (within-compartment Wright-Fisher, fidelity-from-mean-phi coupling, parasite
* advantage, between-compartment selection + random assortment). Observation adds NO RNG draws,
* so trajectories are bit-identical to ribozyme.c. Verified by reproducing its 3-seed table.
*
* MEASURE (following the paper's Materials and Methods):
* Phi_r = I(X_t ; X_t+1) - sum_i I(X_t^i ; X_t+1)
* i.e. Rosas/Mediano-style causal emergence: predictive information held by the WHOLE about its
* own future, beyond what the parts hold about that future. The paper notes the full multivariate
* estimate is intractable and uses the MINIMUM-INFORMATION BIPARTITION to reduce to two
* components; we do the same (exhaustive over bipartitions here, since our part count is small):
* Phi_r = I(X_t ; X_t+1) - I(A_t ; X_t+1) - I(B_t ; X_t+1)
* with (A,B) the bipartition minimising the normalised cross-cut lagged MI
* [ I(A_t;B_t+1) + I(B_t;A_t+1) ] / min(|A|,|B|).
* Substrate: the pool-wide histogram over hairpin score phi (our analogue of GARD's "molecular
* type counts" -- the compositional state), centred-log-ratio transformed with a 0.5 pseudocount,
* last component dropped for full rank -- exactly the paper's preprocessing. MI is estimated
* Gaussian (log-determinant form) with Ledoit-Wolf shrinkage toward the diagonal, so the
* regularisation is data-driven and not a tuned knob.
*
* PRE-REGISTERED PREDICTIONS (written before the first run; report whatever happens):
* P1 EARLY SEPARATION. Computed on the FIRST 25% of generations only (0-49 of 200), Phi_r is
* HIGHER in the arms that later ignite (csz 4, 8) than in the arms that fail (well-mixed,
* csz 16, 32).
* P2 SPIKES PRECEDE. Igniting runs show punctuated Phi_r excursions (>3 sd above that run's own
* mean) occurring BEFORE the ignition generation (first gen at which pool mean phi exceeds
* its start value + 1.0).
* P3 CONTROL -- the prediction that decides whether Phi_r adds anything. Raw mean phi over the
* SAME first 25% does NOT separate the arms. If raw phi already separates them, Phi_r is
* redundant and P1 is worthless even if it passes. This is the analogue of the paper's
* "beat the raw-composition baseline" test, and it is the one that can kill the result.
*
* Any of these may fail. A null is a result. Do not tune to force a positive.
*
* build: cc -O2 -Wall -std=c11 -o ribophi ribophi.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>
/* ---------- model constants: verbatim from ribozyme.c ---------- */
#define L 24
#define PHIMAX (L/2) /* max hairpin pairs = 12 */
#define TOTAL 4800
#define GEN 200
#define F0 0.85
#define F1 0.999
#define BETA 2.0
/* ---------- probe constants ---------- */
#define NBFULL 13 /* phi in 0..12 -> full histogram */
#define NB 8 /* merged bins: <=3,4,5,6,7,8,9,>=10 (primary substrate) */
#define P (NB-1) /* parts after dropping one CLR component */
#define PFULL (NBFULL-1)
#define WIN 50 /* window length in generations = first 25% of GEN */
#define STRIDE 2
#define NSEED 20
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 int phi(uint32_t s){ int c=0; for(int i=0;i<PHIMAX;i++){ int a=(s>>i)&1, b=(s>>(L-1-i))&1; if(a!=b) c++; } return c; }
static uint32_t replicate(uint32_t t,double f){ uint32_t c=0; for(int i=0;i<L;i++){ int b=(t>>i)&1; int nb=(urand()<f)?(1-b):b; if(nb) c|=(1u<<i);} return c; }
static uint32_t pool[TOTAL], nxt[TOTAL];
/* ---------- the model, verbatim + observation (no RNG draws added) ---------- */
static double run(int ncomp,int csz,int group,double *f_end,double *phi_start,
double histf[GEN][NBFULL], double mphi_series[GEN])
{
for(int i=0;i<TOTAL;i++){ uint32_t s=0; for(int b=0;b<L;b++) if(urand()<0.5) s|=(1u<<b); pool[i]=s; }
{ double sp=0; for(int i=0;i<TOTAL;i++) sp+=phi(pool[i]); *phi_start=sp/TOTAL; }
for(int g=0;g<GEN;g++){
for(int cmp=0;cmp<ncomp;cmp++){
int base=cmp*csz; double sphi=0; for(int k=0;k<csz;k++) sphi+=phi(pool[base+k]); double mphi=sphi/csz;
double f=F0+(F1-F0)*(mphi/PHIMAX);
double tot=0; static double w[TOTAL]; static double cum[TOTAL];
for(int k=0;k<csz;k++){ double ww=1.0+BETA*(1.0-(double)phi(pool[base+k])/PHIMAX); tot+=ww; cum[base+k]=tot; (void)w; }
for(int k=0;k<csz;k++){ double r=urand()*tot; int lo=0,hi=csz-1; while(lo<hi){int mid=(lo+hi)/2; if(cum[base+mid]<r)lo=mid+1; else hi=mid;} nxt[base+k]=replicate(pool[base+lo],f); }
}
for(int i=0;i<TOTAL;i++) pool[i]=nxt[i];
if(group){
static double Pv[TOTAL/1], cumP[TOTAL/1]; double tp=0;
for(int cmp=0;cmp<ncomp;cmp++){ int base=cmp*csz; double sphi=0; for(int k=0;k<csz;k++) sphi+=phi(pool[base+k]); tp+=sphi/csz; cumP[cmp]=tp; Pv[cmp]=0; (void)Pv; }
for(int cmp=0;cmp<ncomp;cmp++){
double r=urand()*tp; int lo=0,hi=ncomp-1; while(lo<hi){int mid=(lo+hi)/2; if(cumP[mid]<r)lo=mid+1; else hi=mid;} int par=lo;
for(int k=0;k<csz;k++){ int src=par*csz+(int)(xs()%csz); nxt[cmp*csz+k]=pool[src]; }
}
for(int i=0;i<TOTAL;i++) pool[i]=nxt[i];
}
/* --- observation only: pool-wide phi histogram + mean phi. No RNG consumed. --- */
if(histf){
for(int b=0;b<NBFULL;b++) histf[g][b]=0.0;
double sp=0;
for(int i=0;i<TOTAL;i++){ int p=phi(pool[i]); histf[g][p]+=1.0; sp+=p; }
if(mphi_series) mphi_series[g]=sp/TOTAL;
}
}
double sp=0; for(int i=0;i<TOTAL;i++) sp+=phi(pool[i]); double mphi=sp/TOTAL;
*f_end=F0+(F1-F0)*(mphi/PHIMAX); return mphi;
}
/* ---------- compositional preprocessing: merge bins, CLR, drop last ---------- */
static void merge_bins(const double *full, double *out) /* NBFULL -> NB */
{
for(int b=0;b<NB;b++) out[b]=0.0;
for(int p=0;p<NBFULL;p++){
int b;
if(p<=3) b=0; else if(p>=10) b=NB-1; else b=p-3; /* 4,5,6,7,8,9 -> 1..6 */
out[b]+=full[p];
}
}
static void clr(const double *cnt,int nb,double *out) /* -> nb-1 components */
{
double tot=0; for(int i=0;i<nb;i++) tot+=cnt[i]+0.5;
double lg[NBFULL]; double mean=0;
for(int i=0;i<nb;i++){ lg[i]=log((cnt[i]+0.5)/tot); mean+=lg[i]; }
mean/=nb;
for(int i=0;i<nb-1;i++) out[i]=lg[i]-mean;
}
/* ---------- Gaussian MI via log-determinants, Ledoit-Wolf shrinkage ---------- */
/* Z: n x d row-major (already the joint [past | future] vector). Builds shrunk covariance. */
static void covariance_lw(const double *Z,int n,int d,double *S)
{
double *mu=calloc(d,sizeof(double));
for(int t=0;t<n;t++) for(int j=0;j<d;j++) mu[j]+=Z[t*d+j];
for(int j=0;j<d;j++) mu[j]/=n;
double *Y=malloc((size_t)n*d*sizeof(double));
for(int t=0;t<n;t++) for(int j=0;j<d;j++) Y[t*d+j]=Z[t*d+j]-mu[j];
for(int i=0;i<d;i++) for(int j=0;j<d;j++){
double s=0; for(int t=0;t<n;t++) s+=Y[t*d+i]*Y[t*d+j];
S[i*d+j]=s/n;
}
/* Ledoit-Wolf shrinkage intensity toward the diagonal target */
double num=0, den=0;
for(int i=0;i<d;i++) for(int j=0;j<d;j++){
if(i==j) continue;
double sij=S[i*d+j], v=0;
for(int t=0;t<n;t++){ double e=Y[t*d+i]*Y[t*d+j]-sij; v+=e*e; }
v*= (double)n/((double)(n-1)*(n-1)*(n-1));
num+=v; den+=sij*sij;
}
double lam = (den>0)? num/den : 1.0;
if(lam<0) lam=0;
if(lam>1) lam=1;
for(int i=0;i<d;i++) for(int j=0;j<d;j++) if(i!=j) S[i*d+j]*=(1.0-lam);
/* unbiased-ish rescale to sample covariance convention */
double sc=(double)n/(double)(n-1);
for(int i=0;i<d;i++) for(int j=0;j<d;j++) S[i*d+j]*=sc;
free(mu); free(Y);
}
/* log|S[idx,idx]| via Cholesky, with jitter fallback. Returns 0 and sets *ok=0 on failure. */
static double logdet_sub(const double *S,int d,const int *idx,int k,int *ok)
{
double *M=malloc((size_t)k*k*sizeof(double));
for(int a=0;a<k;a++) for(int b=0;b<k;b++) M[a*k+b]=S[idx[a]*d+idx[b]];
double tr=0; for(int a=0;a<k;a++) tr+=M[a*k+a];
double jitter=1e-10*(tr/k+1e-12);
for(int attempt=0;attempt<8;attempt++){
double *C=malloc((size_t)k*k*sizeof(double));
memcpy(C,M,(size_t)k*k*sizeof(double));
for(int a=0;a<k;a++) C[a*k+a]+=jitter;
int bad=0; double ld=0;
for(int a=0;a<k && !bad;a++){
for(int b=0;b<=a;b++){
double s=C[a*k+b];
for(int m=0;m<b;m++) s-=C[a*k+m]*C[b*k+m];
if(a==b){ if(s<=0){ bad=1; break; } C[a*k+a]=sqrt(s); ld+=2.0*log(C[a*k+a]); }
else C[a*k+b]=s/C[b*k+b];
}
}
free(C);
if(!bad){ free(M); *ok=1; return ld; }
jitter = (jitter>0)? jitter*100.0 : 1e-9*(tr/k+1e-12);
}
free(M); *ok=0; return 0.0;
}
static double mi_sub(const double *S,int d,const int *A,int ka,const int *B,int kb,int *ok)
{
int *U=malloc((size_t)(ka+kb)*sizeof(int));
memcpy(U,A,(size_t)ka*sizeof(int)); memcpy(U+ka,B,(size_t)kb*sizeof(int));
int o1,o2,o3;
double la=logdet_sub(S,d,A,ka,&o1), lb=logdet_sub(S,d,B,kb,&o2), lu=logdet_sub(S,d,U,ka+kb,&o3);
free(U);
*ok = o1&&o2&&o3;
return 0.5*(la+lb-lu);
}
/* Phi_r over n transition pairs of p-dim CLR vectors.
pairs: n x (2p) row-major, [past(0..p-1) | future(p..2p-1)]. Returns 0 and *ok=0 on failure. */
static double phi_r(const double *pairs,int n,int p,int *ok)
{
int d=2*p;
double *S=malloc((size_t)d*d*sizeof(double));
covariance_lw(pairs,n,d,S);
int *past=malloc(p*sizeof(int)), *fut=malloc(p*sizeof(int));
for(int i=0;i<p;i++){ past[i]=i; fut[i]=p+i; }
int o;
double I_whole = mi_sub(S,d,past,p,fut,p,&o);
if(!o){ free(S);free(past);free(fut); *ok=0; return 0.0; }
/* minimum-information bipartition over the p parts */
int bestA[P+1], bestB[P+1], bka=0, bkb=0; double bestcut=1e300; int found=0;
int *Ap=malloc(p*sizeof(int)), *Bp=malloc(p*sizeof(int));
int *Af=malloc(p*sizeof(int)), *Bf=malloc(p*sizeof(int));
for(unsigned m=1; m < (1u<<(p-1)); m++){ /* fix part 0 in A -> unordered bipartitions */
int ka=0,kb=0;
for(int i=0;i<p;i++){
int inA = (i==0)? 1 : ((m>>(i-1))&1);
if(inA){ Ap[ka]=i; Af[ka]=p+i; ka++; } else { Bp[kb]=i; Bf[kb]=p+i; kb++; }
}
if(ka==0||kb==0) continue;
int oa,ob;
double c1=mi_sub(S,d,Ap,ka,Bf,kb,&oa);
double c2=mi_sub(S,d,Bp,kb,Af,ka,&ob);
if(!oa||!ob) continue;
int mn = (ka<kb)?ka:kb;
double cut=(c1+c2)/(double)mn;
if(cut<bestcut){ bestcut=cut; bka=ka; bkb=kb;
memcpy(bestA,Ap,(size_t)ka*sizeof(int)); memcpy(bestB,Bp,(size_t)kb*sizeof(int)); found=1; }
}
free(Ap);free(Bp);free(Af);free(Bf);
if(!found){ free(S);free(past);free(fut); *ok=0; return 0.0; }
int o1,o2;
double I_A = mi_sub(S,d,bestA,bka,fut,p,&o1);
double I_B = mi_sub(S,d,bestB,bkb,fut,p,&o2);
free(S);free(past);free(fut);
*ok = o1&&o2;
return I_whole - I_A - I_B;
}
/* ---------- helpers ---------- */
static int cmpd(const void*a,const void*b){ double x=*(const double*)a,y=*(const double*)b; return (x>y)-(x<y); }
static double median(double *v,int n){ double *c=malloc(n*sizeof(double)); memcpy(c,v,n*sizeof(double)); qsort(c,n,sizeof(double),cmpd); double m=(n%2)?c[n/2]:0.5*(c[n/2-1]+c[n/2]); free(c); return m; }
static double meanv(const double *v,int n){ double s=0; for(int i=0;i<n;i++) s+=v[i]; return s/n; }
static double sdv(const double *v,int n){ double m=meanv(v,n),s=0; for(int i=0;i<n;i++) s+=(v[i]-m)*(v[i]-m); return sqrt(s/(n>1?n-1:1)); }
/* Mann-Whitney U -> normal approximation two-sided p (no ties correction; n=20 per group) */
static double mannwhitney_p(const double *x,int nx,const double *y,int ny)
{
double U=0;
for(int i=0;i<nx;i++) for(int j=0;j<ny;j++){ if(x[i]>y[j]) U+=1.0; else if(x[i]==y[j]) U+=0.5; }
double mu=0.5*nx*ny, sd=sqrt((double)nx*ny*(nx+ny+1)/12.0);
if(sd<=0) return 1.0;
double z=fabs(U-mu)/sd;
return erfc(z/sqrt(2.0));
}
typedef struct { const char*name; int ncomp,csz,group; } Arm;
/* Permutation null: destroy the past->future pairing by shuffling the future block across rows,
recompute Phi_r R times, return mean and sd of the null distribution. */
static void phi_r_null(const double *pairs,int n,int p,int R,double *nmean,double *nsd)
{
double *buf=malloc((size_t)n*2*p*sizeof(double));
double *vals=malloc((size_t)R*sizeof(double));
int *perm=malloc((size_t)n*sizeof(int));
int nv=0;
for(int r=0;r<R;r++){
memcpy(buf,pairs,(size_t)n*2*p*sizeof(double));
for(int i=0;i<n;i++) perm[i]=i;
for(int i=n-1;i>0;i--){ int j=(int)(xs()%(unsigned long long)(i+1)); int tmp=perm[i]; perm[i]=perm[j]; perm[j]=tmp; }
for(int i=0;i<n;i++) for(int k=0;k<p;k++) buf[i*2*p+p+k]=pairs[perm[i]*2*p+p+k];
int ok; double v=phi_r(buf,n,p,&ok);
if(ok) vals[nv++]=v;
}
*nmean = nv? meanv(vals,nv):NAN;
*nsd = nv>1? sdv(vals,nv):NAN;
free(buf);free(vals);free(perm);
}
int main(void)
{
static Arm arms[]={
{"well-mixed", 1, TOTAL, 0},
{"comp 1200x4", 1200, 4, 1},
{"comp 600x8", 600, 8, 1},
{"comp 300x16", 300, 16, 1},
{"comp 150x32", 150, 32, 1},
};
const int NA=(int)(sizeof(arms)/sizeof(arms[0]));
static double histf[GEN][NBFULL];
static double mser[GEN];
/* ============ 0a. estimator validation on synthetic data with KNOWN answers ============
A null result downstream is uninterpretable unless the estimator is shown to respond
with the right SIGN and SENSITIVITY first. Three systems, same p and n as the real probe:
(a) independent -> Phi_r ~ 0
(b) common driver -> Phi_r < 0 (redundancy dominates: each half already predicts)
(c) difference -> Phi_r > 0 (synergy: neither half alone reveals d1-d2) */
printf("=== 0a. VALIDATE estimator on synthetic systems with known answers (p=%d, n=100) ===\n",P);
{
RNG=0xabcdef0123456789ULL;
int n=100; static double sy[100][2*P];
double gz[4];
/* (a) independent */
for(int t=0;t<n;t++) for(int i=0;i<2*P;i++){ double u=urand(),v=urand(); if(u<1e-300)u=1e-300; sy[t][i]=sqrt(-2*log(u))*cos(6.283185307*v); }
int ok; double va=phi_r(&sy[0][0],n,P,&ok);
/* (b) common driver: every part = z + small noise, at both times */
for(int t=0;t<n;t++){
double u=urand(),v=urand(); if(u<1e-300)u=1e-300; double z=sqrt(-2*log(u))*cos(6.283185307*v);
for(int i=0;i<2*P;i++){ double u2=urand(),v2=urand(); if(u2<1e-300)u2=1e-300; double e=sqrt(-2*log(u2))*cos(6.283185307*v2); sy[t][i]=z+0.15*e; }
}
int okb; double vb=phi_r(&sy[0][0],n,P,&okb);
/* (c) synergy under EVERY bipartition. Each part is z + 0.1*d_i, so z swamps it.
future_0 = sum_i w_i d_i with weights w summing to ZERO, so the contrast is invariant
to the common offset z: the WHOLE recovers it exactly (z cancels), while any proper
subset has sum(w) != 0 in general and cannot. Unlike a single synergistic pair, this
cannot be defeated by the MIB choosing a cut that keeps the pair together. */
{
static const double w[P]={+1,-1,+1,-1,+1,-1,0};
double dd[P];
for(int t=0;t<n;t++){
double u=urand(),v=urand(); if(u<1e-300)u=1e-300; double z=sqrt(-2*log(u))*cos(6.283185307*v);
for(int k=0;k<P;k++){ double u2=urand(),v2=urand(); if(u2<1e-300)u2=1e-300; dd[k]=sqrt(-2*log(u2))*cos(6.283185307*v2); }
for(int i=0;i<P;i++) sy[t][i] = z + 0.1*dd[i];
for(int i=0;i<P;i++){ double u3=urand(),v3=urand(); if(u3<1e-300)u3=1e-300; double e=sqrt(-2*log(u3))*cos(6.283185307*v3); sy[t][P+i]=0.05*e; }
double acc=0; for(int i=0;i<P;i++) acc+=w[i]*dd[i];
sy[t][P+0]=acc;
}
}
int okc; double vc=phi_r(&sy[0][0],n,P,&okc);
(void)gz;
/* judge each against its OWN permutation null rather than an arbitrary absolute cut */
double nma,nsa,nmb,nsb,nmc,nsc;
phi_r_null(&sy[0][0],n,P,40,&nmc,&nsc); /* sy currently holds (c) */
double zc = (nsc>0)? (vc-nmc)/nsc : NAN;
printf(" (c) synergy (zero-sum contrast) -> Phi_r = %+8.5f vs shuffled null %+.5f +- %.5f z=%+.1f\n", vc,nmc,nsc,zc);
/* regenerate (a) and (b) to get their nulls too */
for(int t=0;t<n;t++) for(int i=0;i<2*P;i++){ double u=urand(),v2=urand(); if(u<1e-300)u=1e-300; sy[t][i]=sqrt(-2*log(u))*cos(6.283185307*v2); }
phi_r_null(&sy[0][0],n,P,40,&nma,&nsa);
double za=(nsa>0)?(va-nma)/nsa:NAN;
for(int t=0;t<n;t++){
double u=urand(),v2=urand(); if(u<1e-300)u=1e-300; double z=sqrt(-2*log(u))*cos(6.283185307*v2);
for(int i=0;i<2*P;i++){ double u2=urand(),v3=urand(); if(u2<1e-300)u2=1e-300; double e=sqrt(-2*log(u2))*cos(6.283185307*v3); sy[t][i]=z+0.15*e; }
}
phi_r_null(&sy[0][0],n,P,40,&nmb,&nsb);
double zb=(nsb>0)?(vb-nmb)/nsb:NAN;
printf(" (a) independent parts -> Phi_r = %+8.5f vs shuffled null %+.5f +- %.5f z=%+.1f expect ~0\n", va,nma,nsa,za);
printf(" (b) common driver -> Phi_r = %+8.5f vs shuffled null %+.5f +- %.5f z=%+.1f expect <0\n", vb,nmb,nsb,zb);
(void)ok;(void)okb;(void)okc;
int pa=(fabs(za)<3.0), pb=(zb<-3.0), pc=(zc>3.0);
printf(" verdict: (a) %s (b) %s (c) %s\n", pa?"PASS":"FAIL", pb?"PASS":"FAIL", pc?"PASS":"FAIL");
if(pa&&pb&&pc) printf(" => estimator recovers all three signs against its own null. A null below is informative.\n");
else printf(" => ESTIMATOR FAILS A SIGN TEST. Nulls below are confounded and must NOT be read as\n claims about Phi_r itself.\n");
printf(" NOTE: |z| for redundancy (b) far exceeds that for synergy (c) -- this Gaussian estimator\n is much more sensitive to redundancy than to synergy. Carry that caveat into the results.\n\n");
}
/* ================= 0. reproduce ribozyme.c on its own 3 seeds ================= */
printf("=== 0. VERIFY: instrumented model reproduces ribozyme.c (same 3 seeds) ===\n");
printf(" arm | phi start -> end | fidelity end (ribozyme.c reference)\n");
unsigned long long vs[]={0x9e3779b97f4a7c15ULL,0xd1b54a32d192ed03ULL,0x2545f4914f6cdd1dULL};
for(int a=0;a<NA;a++){
double S=0,E=0,F=0;
for(int s=0;s<3;s++){ double fe,ps; RNG=vs[s];
E+=run(arms[a].ncomp,arms[a].csz,arms[a].group,&fe,&ps,histf,mser); S+=ps; F+=fe; }
printf(" %-15s | %.2f -> %5.2f | %.3f\n",arms[a].name,S/3,E/3,F/3);
}
/* ================= 1. run NSEED seeds per arm, collect Phi_r ================= */
printf("\n=== 1. Phi_r probe: %d seeds x %d arms, GEN=%d, window=%d gens (first 25%%) ===\n",NSEED,NA,GEN,WIN);
printf(" substrate = pool-wide phi histogram (%d merged bins -> %d CLR parts), Gaussian MI,\n",NB,P);
printf(" Ledoit-Wolf shrinkage, exhaustive minimum-information bipartition.\n\n");
static double earlyPhi[8][NSEED]; /* Phi_r on gens 0..WIN-1 */
static double earlyRawPhi[8][NSEED]; /* CONTROL: raw mean phi, same gens */
static double finalPhi[8][NSEED]; /* outcome: end mean phi */
static int spikeBefore[8], spikeAny[8], ignited[8];
static int ignGen[8][NSEED]; /* generation at which each run ignites (-1 = never) */
static double clrAll[8][NSEED][GEN][P];/* CLR series kept for the pre-ignition test */
static double mserAll[8][NSEED][GEN];
for(int a=0;a<NA;a++){
spikeBefore[a]=spikeAny[a]=ignited[a]=0;
for(int s=0;s<NSEED;s++){
double fe,ps;
RNG = 0x9e3779b97f4a7c15ULL ^ (0x100000001b3ULL*(unsigned long long)(s+1));
for(int k=0;k<3;k++) xs(); /* decorrelate seed streams */
double endphi = run(arms[a].ncomp,arms[a].csz,arms[a].group,&fe,&ps,histf,mser);
finalPhi[a][s]=endphi;
/* CONTROL: raw compositional signal over the same early window */
double rp=0; for(int g=0;g<WIN;g++) rp+=mser[g];
earlyRawPhi[a][s]=rp/WIN;
/* build CLR series once */
static double clrser[GEN][P];
for(int g=0;g<GEN;g++){ double m8[NB]; merge_bins(histf[g],m8); clr(m8,NB,clrser[g]); }
memcpy(clrAll[a][s],clrser,sizeof(clrser));
memcpy(mserAll[a][s],mser,sizeof(mser));
{ int gg=-1; for(int g=0;g<GEN;g++) if(mser[g] > ps+1.0){ gg=g; break; } ignGen[a][s]=gg; }
/* Phi_r on the early window only */
{
int n=WIN-1; static double pairs[GEN][2*P];
for(int t=0;t<n;t++){ for(int i=0;i<P;i++){ pairs[t][i]=clrser[t][i]; pairs[t][P+i]=clrser[t+1][i]; } }
int ok; double v=phi_r(&pairs[0][0],n,P,&ok);
earlyPhi[a][s] = ok? v : NAN;
}
/* full sliding trajectory, for the spike test */
{
int nw=0; static double traj[GEN]; static int tcenter[GEN];
for(int st=0; st+WIN<=GEN; st+=STRIDE){
int n=WIN-1; static double pairs[GEN][2*P];
for(int t=0;t<n;t++){ for(int i=0;i<P;i++){ pairs[t][i]=clrser[st+t][i]; pairs[t][P+i]=clrser[st+t+1][i]; } }
int ok; double v=phi_r(&pairs[0][0],n,P,&ok);
if(ok){ traj[nw]=v; tcenter[nw]=st+WIN/2; nw++; }
}
(void)meanv;
/* ignition generation: first gen where mean phi exceeds start + 1.0 */
int gign=-1; for(int g=0;g<GEN;g++) if(mser[g] > ps+1.0){ gign=g; break; }
if(gign>=0) ignited[a]++;
if(nw>3){
double mu=meanv(traj,nw), sd=sdv(traj,nw);
for(int i=0;i<nw;i++) if(traj[i] > mu+3.0*sd){
spikeAny[a]++;
if(gign>=0 && tcenter[i]<gign) spikeBefore[a]++;
break;
}
}
}
}
}
/* ================= 2. results ================= */
printf(" arm | end phi | P3 CONTROL raw phi | P1 Phi_r (early, gens 0-%d)\n",WIN-1);
printf(" | (gen%3d)| mean+-sd (0-%2d) | median mean+-sd\n",GEN,WIN-1);
printf(" ----------------+---------+--------------------+---------------------------\n");
for(int a=0;a<NA;a++){
double ep[NSEED]; int ne=0;
for(int s=0;s<NSEED;s++) if(!isnan(earlyPhi[a][s])) ep[ne++]=earlyPhi[a][s];
printf(" %-15s | %6.2f | %5.2f +- %4.2f | %7.4f %7.4f +- %6.4f (n=%d)\n",
arms[a].name, meanv(finalPhi[a],NSEED),
meanv(earlyRawPhi[a],NSEED), sdv(earlyRawPhi[a],NSEED),
ne?median(ep,ne):NAN, ne?meanv(ep,ne):NAN, ne?sdv(ep,ne):NAN, ne);
}
/* pool the igniting arms vs the failing arms */
double ig[2*NSEED], fl[3*NSEED]; int nig=0,nfl=0;
double igR[2*NSEED], flR[3*NSEED]; int nigR=0,nflR=0;
for(int a=0;a<NA;a++){
int isIgnite = (arms[a].csz==4 || arms[a].csz==8) && arms[a].group;
for(int s=0;s<NSEED;s++){
if(!isnan(earlyPhi[a][s])){ if(isIgnite) ig[nig++]=earlyPhi[a][s]; else fl[nfl++]=earlyPhi[a][s]; }
if(isIgnite) igR[nigR++]=earlyRawPhi[a][s]; else flR[nflR++]=earlyRawPhi[a][s];
}
}
printf("\n --- P1: early Phi_r, igniting arms (csz 4,8) vs failing arms (well-mixed,16,32) ---\n");
printf(" igniting: mean %.4f (n=%d) failing: mean %.4f (n=%d) Mann-Whitney p=%.4g\n",
meanv(ig,nig),nig, meanv(fl,nfl),nfl, mannwhitney_p(ig,nig,fl,nfl));
printf(" --- P3 CONTROL: raw mean phi over the SAME early window ---\n");
printf(" igniting: mean %.4f (n=%d) failing: mean %.4f (n=%d) Mann-Whitney p=%.4g\n",
meanv(igR,nigR),nigR, meanv(flR,nflR),nflR, mannwhitney_p(igR,nigR,flR,nflR));
printf("\n --- P2: spikes (>3 sd above run mean) and whether they precede ignition ---\n");
printf(" arm | runs ignited | runs with a spike | spike BEFORE ignition\n");
for(int a=0;a<NA;a++)
printf(" %-15s | %5d / %-5d | %6d / %-8d | %d\n",
arms[a].name, ignited[a],NSEED, spikeAny[a],NSEED, spikeBefore[a]);
/* ================= 3. WHEN does ignition happen? ================= */
printf("\n=== 3. Ignition timing (first gen where pool mean phi > start+1.0) ===\n");
printf(" arm | ignited | ignition gen: min / median / max\n");
int gmin_global=GEN;
for(int a=0;a<NA;a++){
double gv[NSEED]; int ng=0; int mn=GEN,mx=-1;
for(int s=0;s<NSEED;s++) if(ignGen[a][s]>=0){ gv[ng++]=ignGen[a][s]; if(ignGen[a][s]<mn)mn=ignGen[a][s]; if(ignGen[a][s]>mx)mx=ignGen[a][s]; }
if(ng){ printf(" %-15s | %2d/%2d | %d / %.1f / %d\n",arms[a].name,ng,NSEED,mn,median(gv,ng),mx);
if(mn<gmin_global) gmin_global=mn; }
else printf(" %-15s | %2d/%2d | -\n",arms[a].name,ng,NSEED);
}
printf(" => earliest ignition anywhere: gen %d. The window used in section 1 was gens 0-%d,\n",gmin_global,WIN-1);
printf(" so that window was NOT pre-transition. Re-testing strictly before gen %d.\n",gmin_global);
/* ================= 4. strictly PRE-IGNITION window ================= */
int PRE = gmin_global; /* gens 0..PRE-1 contain no ignited state in ANY run */
printf("\n=== 4. P1 re-test on a strictly PRE-IGNITION window (gens 0-%d) ===\n",PRE-1);
if(PRE < 4){
printf(" Only %d pre-ignition generations exist -- too few transitions per run for a per-run\n",PRE);
printf(" estimate. Pooling transition pairs across the %d seeds of each arm instead.\n",NSEED);
}
{
int perRun = PRE-1; /* transition pairs available per run */
if(perRun < 1){
printf(" NO pre-ignition transitions exist at generation resolution. Test cannot be run here.\n");
} else {
printf(" arm | pooled n | Phi_r (pre-ignition) | CONTROL raw phi mean+-sd\n");
printf(" ----------------+----------+----------------------+-------------------------\n");
static double poolPairs[NSEED*GEN][2*P];
double preRawIg[NSEED*8]; int nRawIg=0; double preRawFl[NSEED*8]; int nRawFl=0;
double phirIg[8]; double phirFl[8]; int nIgA=0,nFlA=0;
for(int a=0;a<NA;a++){
int n=0;
for(int s=0;s<NSEED;s++)
for(int t=0;t<perRun;t++){
for(int i=0;i<P;i++){ poolPairs[n][i]=clrAll[a][s][t][i]; poolPairs[n][P+i]=clrAll[a][s][t+1][i]; }
n++;
}
int ok; double v=phi_r(&poolPairs[0][0],n,P,&ok);
double nm,ns; phi_r_null(&poolPairs[0][0],n,P,40,&nm,&ns);
double zz=(ns>0)?(v-nm)/ns:NAN;
double rr[NSEED];
for(int s=0;s<NSEED;s++){ double q=0; for(int t=0;t<PRE;t++) q+=mserAll[a][s][t]; rr[s]=q/PRE; }
printf(" %-15s | %8d | %9.5f (z=%+5.1f) | %.3f +- %.3f\n",
arms[a].name,n, ok? v : NAN, zz, meanv(rr,NSEED), sdv(rr,NSEED));
int isIgnite = (arms[a].csz==4 || arms[a].csz==8) && arms[a].group;
if(isIgnite){ if(ok) phirIg[nIgA++]=v; for(int s=0;s<NSEED;s++) preRawIg[nRawIg++]=rr[s]; }
else { if(ok) phirFl[nFlA++]=v; for(int s=0;s<NSEED;s++) preRawFl[nRawFl++]=rr[s]; }
}
printf("\n P1 (pre-ignition): mean Phi_r igniting arms %.5f vs failing arms %.5f\n",
meanv(phirIg,nIgA), meanv(phirFl,nFlA));
printf(" P3 CONTROL (pre-ignition): raw phi igniting %.4f vs failing %.4f Mann-Whitney p=%.4g\n",
meanv(preRawIg,nRawIg), meanv(preRawFl,nRawFl), mannwhitney_p(preRawIg,nRawIg,preRawFl,nRawFl));
}
}
printf("\n=== verdict is read off the pre-registered predictions above; see ribophi.md ===\n");
return 0;
}