/* rung34phi.c — Phi_r (causal emergence) probe on rung34.c's compartmentalization emergence.
*
* WHY THIS MODEL. The ribozyme probe (ribophi.c) returned a null with a diagnosis: ribozyme.c
* ignites at generation 6-10, so it has NO quiescent pre-replicative phase and cannot host the
* paper's "predict the transition from the early trajectory" protocol at all. rung34.c is the
* better test bed: the compartmentalization trait c starts at 0.10 and climbs to ~0.94 over 500
* generations in mutation steps of 0.03, so a genuine lag phase should exist.
*
* WHAT IS COPIED VERBATIM from rung34.c: RNG, Dof(), binom(), randn(), and the whole run() loop
* (within-compartment parasite growth, cooperator->parasite mutation, trait mutation, daughter
* assortment, survivor selection weighted by cooperator FRACTION). Observation adds no RNG draws.
* Verified by reproducing rung34.c's single-seed result.
*
* SUBSTRATE. Histogram over cooperator FRACTION across the NCOMP compartments (8 equal bins on
* [0,1]), CLR-transformed with a 0.5 pseudocount, last component dropped -- the analogue of GARD's
* molecular type counts. Note the trait histogram would be degenerate in the frozen control (c is
* constant by construction), which is why cooperator fraction is the substrate, not c.
*
* PRE-REGISTERED (written before the first run):
* Q1 WITHIN-ARM (this is the paper's actual design: one parameter set, many seeds, predict
* timing). In the c-evolvable arm, Phi_r computed on a strictly PRE-IGNITION window correlates
* across seeds with the generation at which the trait climb ignites.
* Q2 SPIKES. Igniting runs show Phi_r excursions (>3 sd above the run's own mean) before ignition.
* Q3 CONTROL -- the one that can kill it. Raw mean cooperator fraction over the SAME window does
* NOT predict ignition generation. If it does, Phi_r adds nothing.
* Q4 BETWEEN-ARM. Early Phi_r separates c-evolvable from the frozen control. (Weaker test: the
* arms differ mechanically from generation 0, so a difference here is not decisive.)
*
* Any of these may fail. A null is a result. Do not tune to force a positive.
*
* build: cc -O2 -Wall -std=c11 -o rung34phi rung34phi.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "phir.h"
/* ---------- model constants: verbatim from rung34.c ---------- */
#define NCOMP 600
#define GEN 500
#define SADV 0.6
#define DMIN 2
#define DMAX 40
#define MUTC 0.03
#define MUP 0.10
/* ---------- probe constants ---------- */
#define NB 8 /* cooperator-fraction bins */
#define P (NB-1)
#define NSEED 40
#define STRIDE 5
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 double randn(void){ double u=urand(),v=urand(); if(u<1e-300)u=1e-300; return sqrt(-2*log(u))*cos(6.283185307*v); }
static int Dof(double c){ int D=(int)(DMIN+(1.0-c)*(DMAX-DMIN)+0.5); return D<DMIN?DMIN:(D>DMAX?DMAX:D); }
static int binom(int n,double p){ int k=0; for(int i=0;i<n;i++) if(urand()<p) k++; return k; }
static double cc[NCOMP]; static int nC[NCOMP];
static double dcc[2*NCOMP]; static int dnC[2*NCOMP], dD[2*NCOMP]; static double cum[2*NCOMP];
/* verbatim rung34 dynamics + observation */
static double run(int evolve_c,double *coopf,double cser[GEN],double fser[GEN],double hist[GEN][NB],double histc[GEN][NB])
{
for(int i=0;i<NCOMP;i++){ cc[i]=0.10; nC[i]=Dof(0.10)/2; }
for(int g=0;g<GEN;g++){
for(int i=0;i<NCOMP;i++){
int D=Dof(cc[i]); int nc=nC[i]; if(nc>D)nc=D; int np=D-nc;
double wc=nc, wp=np*(1.0+SADV); double p=(wc+wp>0)?wc/(wc+wp):0;
int g2=binom(2*D,p); g2=binom(g2,1.0-MUP); double frac=(double)g2/(2*D);
for(int d=0;d<2;d++){
double c2 = evolve_c ? cc[i]+MUTC*randn() : cc[i];
if(c2<0)c2=0;
if(c2>1)c2=1;
int D2=Dof(c2);
dcc[2*i+d]=c2; dD[2*i+d]=D2; dnC[2*i+d]=binom(D2,frac);
}
}
double tot=0; for(int i=0;i<2*NCOMP;i++){ tot+=(double)dnC[i]/dD[i]; cum[i]=tot; }
if(tot<=0){ for(int i=0;i<NCOMP;i++){ nC[i]=0; } }
else for(int i=0;i<NCOMP;i++){
double r=urand()*tot; int lo=0,hi=2*NCOMP-1;
while(lo<hi){ int mid=(lo+hi)/2; if(cum[mid]<r) lo=mid+1; else hi=mid; }
cc[i]=dcc[lo]; nC[i]=dnC[lo];
}
/* --- observation only, no RNG consumed --- */
{
double mc=0; long sc=0,sd=0;
for(int b=0;b<NB;b++) hist[g][b]=0.0;
for(int i=0;i<NCOMP;i++){
int D=Dof(cc[i]); mc+=cc[i]; sc+=nC[i]; sd+=D;
double fr=(D>0)?(double)nC[i]/D:0.0;
int b=(int)(fr*NB); if(b>=NB)b=NB-1; if(b<0)b=0;
hist[g][b]+=1.0;
}
if(cser) cser[g]=mc/NCOMP;
if(fser) fser[g]=(double)sc/sd;
}
/* second substrate: distribution of the TRAIT c itself (degenerate in the frozen arm,
so usable only for the within-arm test, where it is the variable that actually moves) */
if(histc){
for(int b=0;b<NB;b++) histc[g][b]=0.0;
for(int i=0;i<NCOMP;i++){ int b=(int)(cc[i]*NB); if(b>=NB)b=NB-1; if(b<0)b=0; histc[g][b]+=1.0; }
}
}
double mc=0; long sc=0,sd=0; for(int i=0;i<NCOMP;i++){ mc+=cc[i]; sc+=nC[i]; sd+=Dof(cc[i]); }
*coopf=(double)sc/sd; return mc/NCOMP;
}
static double hist_all[2][NSEED][GEN][NB];
static double cAll[2][NSEED][GEN], fAll[2][NSEED][GEN];
static double clrAll[2][NSEED][GEN][P];
static double histc_all[2][NSEED][GEN][NB];
static double clrC[2][NSEED][GEN][P];
static int ignGen[2][NSEED];
static double finalC[2][NSEED], finalF[2][NSEED];
int main(void)
{
phir_seed(0x51ed270b0f8e1c4dULL);
/* ---------- 0. reproduce rung34.c on its own seed ---------- */
printf("=== 0. VERIFY: instrumented model reproduces rung34.c (its own seed) ===\n");
{
static double cs[GEN],fs[GEN]; static double hh[GEN][NB],hc[GEN][NB]; double cf;
RNG=0x9e3779b97f4a7c15ULL; double c1=run(1,&cf,cs,fs,hh,hc);
printf(" c EVOLVABLE : final mean c = %.3f (D=%d) cooperation frac = %.3f\n",c1,Dof(c1),cf);
RNG=0x9e3779b97f4a7c15ULL; double c0=run(0,&cf,cs,fs,hh,hc);
printf(" c FROZEN : final mean c = %.3f (D=%d) cooperation frac = %.3f\n",c0,Dof(c0),cf);
printf(" (rung34.c reference: evolvable c ~0.94 with cooperation rescued; frozen collapses)\n\n");
}
/* ---------- 1. run both arms, NSEED seeds ---------- */
const char *armname[2]={"c EVOLVABLE","c FROZEN "};
for(int a=0;a<2;a++){
for(int s=0;s<NSEED;s++){
double cf;
RNG = 0x9e3779b97f4a7c15ULL ^ (0x100000001b3ULL*(unsigned long long)(s+1));
for(int k=0;k<3;k++) xs();
finalC[a][s]=run(a==0?1:0,&cf,cAll[a][s],fAll[a][s],hist_all[a][s],histc_all[a][s]);
finalF[a][s]=cf;
for(int g=0;g<GEN;g++){ phir_clr(hist_all[a][s][g],NB,clrAll[a][s][g]); phir_clr(histc_all[a][s][g],NB,clrC[a][s][g]); }
int gg=-1; for(int g=0;g<GEN;g++) if(cAll[a][s][g] > 0.10+0.20){ gg=g; break; }
ignGen[a][s]=gg;
}
}
/* ---------- 2. ignition timing: is there a lag phase at all? ---------- */
printf("=== 1. Ignition timing (first gen where mean trait c > start+0.20) ===\n");
int gmin=GEN;
for(int a=0;a<2;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(" %s | ignited %2d/%2d | gen min/median/max = %d / %.1f / %d\n",armname[a],ng,NSEED,mn,median(gv,ng),mx);
if(mn<gmin) gmin=mn; }
else printf(" %s | ignited %2d/%2d | (never)\n",armname[a],ng,NSEED);
printf(" final c = %.3f +- %.3f final coop frac = %.3f +- %.3f\n",
meanv(finalC[a],NSEED),sdv(finalC[a],NSEED),meanv(finalF[a],NSEED),sdv(finalF[a],NSEED));
}
if(gmin>=GEN){ printf(" no ignition anywhere -- cannot run the test.\n"); return 1; }
int PRE = gmin; /* gens 0..PRE-1 are pre-ignition in EVERY run */
printf(" => earliest ignition anywhere: gen %d. Pre-ignition window = gens 0-%d (%d transitions/run).\n\n",
gmin,PRE-1,PRE-1);
/* ---------- 3. Q1/Q3: within-arm, does early Phi_r predict ignition timing? ---------- */
printf("=== 2. Q1 WITHIN-ARM (c evolvable): does pre-ignition Phi_r predict ignition gen? ===\n");
{
int n=PRE-1;
if(n<8){
printf(" only %d pre-ignition transitions per run -- too few for a per-run estimate.\n",n);
} else {
static double pairs[GEN][2*P];
double ephi[NSEED], eraw[NSEED], gign[NSEED], fc[NSEED]; int m=0;
for(int s=0;s<NSEED;s++){
if(ignGen[0][s]<0) continue;
for(int t=0;t<n;t++) for(int i=0;i<P;i++){ pairs[t][i]=clrAll[0][s][t][i]; pairs[t][P+i]=clrAll[0][s][t+1][i]; }
int ok; double v=phi_r(&pairs[0][0],n,P,&ok);
if(!ok) continue;
double q=0; for(int t=0;t<PRE;t++) q+=fAll[0][s][t];
ephi[m]=v; eraw[m]=q/PRE; gign[m]=ignGen[0][s]; fc[m]=finalC[0][s]; m++;
}
printf(" n=%d runs, window gens 0-%d\n",m,PRE-1);
printf(" pre-ignition Phi_r : %.5f +- %.5f\n",meanv(ephi,m),sdv(ephi,m));
double r1=spearman(ephi,gign,m), r2=spearman(ephi,fc,m);
double c1=spearman(eraw,gign,m), c2=spearman(eraw,fc,m);
printf(" Q1 Phi_r vs ignition gen : rho = %+.3f (p=%.3g)\n",r1,spearman_p(r1,m));
printf(" Phi_r vs final c : rho = %+.3f (p=%.3g)\n",r2,spearman_p(r2,m));
printf(" Q3 raw coopfr vs ignition gen : rho = %+.3f (p=%.3g) <- CONTROL\n",c1,spearman_p(c1,m));
printf(" raw coopfr vs final c : rho = %+.3f (p=%.3g) <- CONTROL\n",c2,spearman_p(c2,m));
printf(" => Phi_r earns its keep only if |rho| for Phi_r exceeds the control's.\n\n");
/* Same test on the substrate that actually carries the transition: the TRAIT
distribution. Cooperator fraction may simply be the wrong variable to watch. */
printf(" --- same test, substrate = TRAIT c distribution (the variable that moves) ---\n");
double ephiC[NSEED], erawC[NSEED], gignC[NSEED]; int mc2=0;
for(int s=0;s<NSEED;s++){
if(ignGen[0][s]<0) continue;
for(int t=0;t<n;t++) for(int i=0;i<P;i++){ pairs[t][i]=clrC[0][s][t][i]; pairs[t][P+i]=clrC[0][s][t+1][i]; }
int ok; double v=phi_r(&pairs[0][0],n,P,&ok);
if(!ok) continue;
double q=0; for(int t=0;t<PRE;t++) q+=cAll[0][s][t];
ephiC[mc2]=v; erawC[mc2]=q/PRE; gignC[mc2]=ignGen[0][s]; mc2++;
}
double rC=spearman(ephiC,gignC,mc2), cC=spearman(erawC,gignC,mc2);
printf(" pre-ignition Phi_r (trait) : %.5f +- %.5f (n=%d)\n",meanv(ephiC,mc2),sdv(ephiC,mc2),mc2);
printf(" Q1' Phi_r(trait) vs ignition gen : rho = %+.3f (p=%.3g)\n",rC,spearman_p(rC,mc2));
printf(" Q3' raw mean c vs ignition gen : rho = %+.3f (p=%.3g) <- CONTROL\n\n",cC,spearman_p(cC,mc2));
}
}
/* ---------- 4. Q4: between-arm on the pre-ignition window ---------- */
printf("=== 3. Q4 BETWEEN-ARM on the pre-ignition window (gens 0-%d) ===\n",PRE-1);
{
int n=PRE-1;
if(n>=8){
static double pairs[GEN][2*P];
double ev[NSEED],fz[NSEED]; int ne=0,nf=0;
double evR[NSEED],fzR[NSEED];
for(int a=0;a<2;a++) for(int s=0;s<NSEED;s++){
for(int t=0;t<n;t++) for(int i=0;i<P;i++){ pairs[t][i]=clrAll[a][s][t][i]; pairs[t][P+i]=clrAll[a][s][t+1][i]; }
int ok; double v=phi_r(&pairs[0][0],n,P,&ok);
double q=0; for(int t=0;t<PRE;t++) q+=fAll[a][s][t]; q/=PRE;
if(!ok) continue;
if(a==0){ ev[ne]=v; evR[ne]=q; ne++; } else { fz[nf]=v; fzR[nf]=q; nf++; }
}
printf(" Phi_r evolvable %.5f +- %.5f | frozen %.5f +- %.5f | MW p=%.4g\n",
meanv(ev,ne),sdv(ev,ne),meanv(fz,nf),sdv(fz,nf),mannwhitney_p(ev,ne,fz,nf));
printf(" raw coopfr evolvable %.5f +- %.5f | frozen %.5f +- %.5f | MW p=%.4g <- CONTROL\n",
meanv(evR,ne),sdv(evR,ne),meanv(fzR,nf),sdv(fzR,nf),mannwhitney_p(evR,ne,fzR,nf));
}
}
/* ---------- 5. Q2: do spikes precede ignition? ---------- */
printf("\n=== 4. Q2 SPIKES: Phi_r excursions >3 sd above run mean, before ignition ===\n");
{
int WIN=PRE>25?PRE:25;
printf(" sliding window = %d gens, stride %d\n",WIN,STRIDE);
for(int a=0;a<2;a++){
int any=0,before=0,nr=0;
for(int s=0;s<NSEED;s++){
static double pairs[GEN][2*P]; static double traj[GEN]; static int tc[GEN];
int nw=0;
for(int st=0; st+WIN<=GEN; st+=STRIDE){
int n=WIN-1;
for(int t=0;t<n;t++) for(int i=0;i<P;i++){ pairs[t][i]=clrAll[a][s][st+t][i]; pairs[t][P+i]=clrAll[a][s][st+t+1][i]; }
int ok; double v=phi_r(&pairs[0][0],n,P,&ok);
if(ok){ traj[nw]=v; tc[nw]=st+WIN/2; nw++; }
}
if(nw<4) continue;
nr++;
double mu=meanv(traj,nw), sd=sdv(traj,nw);
for(int i=0;i<nw;i++) if(traj[i]>mu+3.0*sd){
any++;
if(ignGen[a][s]>=0 && tc[i]<ignGen[a][s]) before++;
break;
}
}
printf(" %s | runs %2d | with a spike %2d | spike BEFORE ignition %2d\n",armname[a],nr,any,before);
}
}
printf("\n=== verdict read off Q1-Q4 above ===\n");
return 0;
}