| … | … | @@ -35,3 +35,5 @@ coemerge |
| 35 | 35 | integrated |
| 36 | 36 | belief_field_agent |
| 37 | 37 | agent_learn |
| 38 | ribophi | |
| 39 | rung34phi |
| … | … | @@ -0,0 +1,211 @@ |
| 1 | /* phir.h — Phi_r (PhiID causal emergence) estimator, shared by the probe programs. | |
| 2 | * | |
| 3 | * Phi_r = I(X_t ; X_t+1) - sum_i I(X_t^i ; X_t+1) | |
| 4 | * i.e. predictive information the WHOLE holds about its own future, beyond what its parts hold | |
| 5 | * about that future. Following Pigozzi & Levin (arXiv 2607.28250) Materials and Methods, the | |
| 6 | * intractable multivariate form is reduced via the MINIMUM-INFORMATION BIPARTITION to | |
| 7 | * Phi_r = I(X_t ; X_t+1) - I(A_t ; X_t+1) - I(B_t ; X_t+1) | |
| 8 | * with (A,B) minimising the normalised cross-cut lagged MI [I(A_t;B_t+1)+I(B_t;A_t+1)]/min(|A|,|B|). | |
| 9 | * | |
| 10 | * MI is estimated Gaussian (log-determinant form) with Ledoit-Wolf shrinkage toward the diagonal, | |
| 11 | * so the regularisation is data-driven rather than a tuned knob. | |
| 12 | * | |
| 13 | * KNOWN LIMITATION, measured: this estimator is far more sensitive to redundancy than to synergy | |
| 14 | * (|z| ~ 165 vs ~3.4 on matched synthetic controls). Carry that caveat into any result. | |
| 15 | * | |
| 16 | * The permutation null uses its own RNG so that null draws never perturb a caller's model stream. | |
| 17 | */ | |
| 18 | #ifndef PHIR_H | |
| 19 | #define PHIR_H | |
| 20 | ||
| 21 | #include <math.h> | |
| 22 | #include <stdlib.h> | |
| 23 | #include <string.h> | |
| 24 | ||
| 25 | #define PHIR_MAXP 24 | |
| 26 | ||
| 27 | static unsigned long long PHIR_RNG = 0x243f6a8885a308d3ULL; | |
| 28 | static unsigned long long phir_xs(void){ PHIR_RNG^=PHIR_RNG<<13; PHIR_RNG^=PHIR_RNG>>7; PHIR_RNG^=PHIR_RNG<<17; return PHIR_RNG; } | |
| 29 | static void phir_seed(unsigned long long s){ PHIR_RNG = s? s : 0x243f6a8885a308d3ULL; } | |
| 30 | ||
| 31 | /* ---------- generic stats ---------- */ | |
| 32 | static int phir_cmpd(const void*a,const void*b){ double x=*(const double*)a,y=*(const double*)b; return (x>y)-(x<y); } | |
| 33 | static double meanv(const double *v,int n){ double s=0; for(int i=0;i<n;i++) s+=v[i]; return n? s/n : NAN; } | |
| 34 | 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)); } | |
| 35 | static double median(double *v,int n){ double *c=malloc((size_t)n*sizeof(double)); memcpy(c,v,(size_t)n*sizeof(double)); qsort(c,n,sizeof(double),phir_cmpd); double m=(n%2)?c[n/2]:0.5*(c[n/2-1]+c[n/2]); free(c); return m; } | |
| 36 | ||
| 37 | /* Mann-Whitney U -> two-sided normal-approximation p (no ties correction) */ | |
| 38 | static double mannwhitney_p(const double *x,int nx,const double *y,int ny) | |
| 39 | { | |
| 40 | double U=0; | |
| 41 | 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; } | |
| 42 | double mu=0.5*nx*ny, sd=sqrt((double)nx*ny*(nx+ny+1)/12.0); | |
| 43 | if(sd<=0) return 1.0; | |
| 44 | return erfc(fabs(U-mu)/sd/sqrt(2.0)); | |
| 45 | } | |
| 46 | /* Spearman rank correlation (average ranks for ties) */ | |
| 47 | static void phir_ranks(const double *v,int n,double *r) | |
| 48 | { | |
| 49 | int *idx=malloc((size_t)n*sizeof(int)); | |
| 50 | for(int i=0;i<n;i++) idx[i]=i; | |
| 51 | for(int i=1;i<n;i++){ int k=idx[i],j=i-1; while(j>=0 && v[idx[j]]>v[k]){ idx[j+1]=idx[j]; j--; } idx[j+1]=k; } | |
| 52 | int i=0; | |
| 53 | while(i<n){ int j=i; while(j+1<n && v[idx[j+1]]==v[idx[i]]) j++; | |
| 54 | double avg=0.5*((double)i+(double)j)+1.0; | |
| 55 | for(int k=i;k<=j;k++) r[idx[k]]=avg; | |
| 56 | i=j+1; } | |
| 57 | free(idx); | |
| 58 | } | |
| 59 | static double spearman(const double *x,const double *y,int n) | |
| 60 | { | |
| 61 | double *rx=malloc((size_t)n*sizeof(double)), *ry=malloc((size_t)n*sizeof(double)); | |
| 62 | phir_ranks(x,n,rx); phir_ranks(y,n,ry); | |
| 63 | double mx=meanv(rx,n),my=meanv(ry,n),sxy=0,sxx=0,syy=0; | |
| 64 | for(int i=0;i<n;i++){ double a=rx[i]-mx,b=ry[i]-my; sxy+=a*b; sxx+=a*a; syy+=b*b; } | |
| 65 | free(rx);free(ry); | |
| 66 | return (sxx>0&&syy>0)? sxy/sqrt(sxx*syy) : NAN; | |
| 67 | } | |
| 68 | static double spearman_p(double rho,int n) /* two-sided, t approximation */ | |
| 69 | { | |
| 70 | if(n<4||isnan(rho)) return NAN; | |
| 71 | double r=rho; if(r>0.999999)r=0.999999; if(r<-0.999999)r=-0.999999; | |
| 72 | double t=fabs(r)*sqrt((n-2)/(1-r*r)); | |
| 73 | return erfc(t/sqrt(2.0)); /* normal approx; conservative enough for n>=20 */ | |
| 74 | } | |
| 75 | ||
| 76 | /* ---------- compositional preprocessing: CLR with pseudocount, drop last component ---------- */ | |
| 77 | static void phir_clr(const double *cnt,int nb,double *out) | |
| 78 | { | |
| 79 | double tot=0; for(int i=0;i<nb;i++) tot+=cnt[i]+0.5; | |
| 80 | double *lg=malloc((size_t)nb*sizeof(double)); double mean=0; | |
| 81 | for(int i=0;i<nb;i++){ lg[i]=log((cnt[i]+0.5)/tot); mean+=lg[i]; } | |
| 82 | mean/=nb; | |
| 83 | for(int i=0;i<nb-1;i++) out[i]=lg[i]-mean; | |
| 84 | free(lg); | |
| 85 | } | |
| 86 | ||
| 87 | /* ---------- Gaussian MI via log-determinants, Ledoit-Wolf shrinkage ---------- */ | |
| 88 | static void phir_cov_lw(const double *Z,int n,int d,double *S) | |
| 89 | { | |
| 90 | double *mu=calloc((size_t)d,sizeof(double)); | |
| 91 | for(int t=0;t<n;t++) for(int j=0;j<d;j++) mu[j]+=Z[(size_t)t*d+j]; | |
| 92 | for(int j=0;j<d;j++) mu[j]/=n; | |
| 93 | double *Y=malloc((size_t)n*d*sizeof(double)); | |
| 94 | for(int t=0;t<n;t++) for(int j=0;j<d;j++) Y[(size_t)t*d+j]=Z[(size_t)t*d+j]-mu[j]; | |
| 95 | for(int i=0;i<d;i++) for(int j=0;j<d;j++){ | |
| 96 | double s=0; for(int t=0;t<n;t++) s+=Y[(size_t)t*d+i]*Y[(size_t)t*d+j]; | |
| 97 | S[i*d+j]=s/n; | |
| 98 | } | |
| 99 | double num=0, den=0; | |
| 100 | for(int i=0;i<d;i++) for(int j=0;j<d;j++){ | |
| 101 | if(i==j) continue; | |
| 102 | double sij=S[i*d+j], v=0; | |
| 103 | for(int t=0;t<n;t++){ double e=Y[(size_t)t*d+i]*Y[(size_t)t*d+j]-sij; v+=e*e; } | |
| 104 | v *= (double)n/((double)(n-1)*(n-1)*(n-1)); | |
| 105 | num+=v; den+=sij*sij; | |
| 106 | } | |
| 107 | double lam = (den>0)? num/den : 1.0; | |
| 108 | if(lam<0) lam=0; | |
| 109 | if(lam>1) lam=1; | |
| 110 | for(int i=0;i<d;i++) for(int j=0;j<d;j++) if(i!=j) S[i*d+j]*=(1.0-lam); | |
| 111 | double sc=(double)n/(double)(n-1); | |
| 112 | for(int i=0;i<d;i++) for(int j=0;j<d;j++) S[i*d+j]*=sc; | |
| 113 | free(mu); free(Y); | |
| 114 | } | |
| 115 | static double phir_logdet(const double *S,int d,const int *idx,int k,int *ok) | |
| 116 | { | |
| 117 | double *M=malloc((size_t)k*k*sizeof(double)); | |
| 118 | for(int a=0;a<k;a++) for(int b=0;b<k;b++) M[a*k+b]=S[idx[a]*d+idx[b]]; | |
| 119 | double tr=0; for(int a=0;a<k;a++) tr+=M[a*k+a]; | |
| 120 | double jitter=1e-10*(tr/k+1e-12); | |
| 121 | for(int attempt=0;attempt<8;attempt++){ | |
| 122 | double *C=malloc((size_t)k*k*sizeof(double)); | |
| 123 | memcpy(C,M,(size_t)k*k*sizeof(double)); | |
| 124 | for(int a=0;a<k;a++) C[a*k+a]+=jitter; | |
| 125 | int bad=0; double ld=0; | |
| 126 | for(int a=0;a<k && !bad;a++){ | |
| 127 | for(int b=0;b<=a;b++){ | |
| 128 | double s=C[a*k+b]; | |
| 129 | for(int m=0;m<b;m++) s-=C[a*k+m]*C[b*k+m]; | |
| 130 | if(a==b){ if(s<=0){ bad=1; break; } C[a*k+a]=sqrt(s); ld+=2.0*log(C[a*k+a]); } | |
| 131 | else C[a*k+b]=s/C[b*k+b]; | |
| 132 | } | |
| 133 | } | |
| 134 | free(C); | |
| 135 | if(!bad){ free(M); *ok=1; return ld; } | |
| 136 | jitter = (jitter>0)? jitter*100.0 : 1e-9*(tr/k+1e-12); | |
| 137 | } | |
| 138 | free(M); *ok=0; return 0.0; | |
| 139 | } | |
| 140 | static double phir_mi(const double *S,int d,const int *A,int ka,const int *B,int kb,int *ok) | |
| 141 | { | |
| 142 | int *U=malloc((size_t)(ka+kb)*sizeof(int)); | |
| 143 | memcpy(U,A,(size_t)ka*sizeof(int)); memcpy(U+ka,B,(size_t)kb*sizeof(int)); | |
| 144 | int o1,o2,o3; | |
| 145 | double la=phir_logdet(S,d,A,ka,&o1), lb=phir_logdet(S,d,B,kb,&o2), lu=phir_logdet(S,d,U,ka+kb,&o3); | |
| 146 | free(U); | |
| 147 | *ok=o1&&o2&&o3; | |
| 148 | return 0.5*(la+lb-lu); | |
| 149 | } | |
| 150 | ||
| 151 | /* pairs: n x (2p) row-major, [past(0..p-1) | future(p..2p-1)] */ | |
| 152 | static double phi_r(const double *pairs,int n,int p,int *ok) | |
| 153 | { | |
| 154 | int d=2*p; | |
| 155 | double *S=malloc((size_t)d*d*sizeof(double)); | |
| 156 | phir_cov_lw(pairs,n,d,S); | |
| 157 | int past[PHIR_MAXP], fut[PHIR_MAXP]; | |
| 158 | for(int i=0;i<p;i++){ past[i]=i; fut[i]=p+i; } | |
| 159 | int o; | |
| 160 | double I_whole=phir_mi(S,d,past,p,fut,p,&o); | |
| 161 | if(!o){ free(S); *ok=0; return 0.0; } | |
| 162 | ||
| 163 | int bestA[PHIR_MAXP], bestB[PHIR_MAXP], bka=0, bkb=0, found=0; | |
| 164 | double bestcut=1e300; | |
| 165 | int Ap[PHIR_MAXP],Bp[PHIR_MAXP],Af[PHIR_MAXP],Bf[PHIR_MAXP]; | |
| 166 | for(unsigned m=1; m<(1u<<(p-1)); m++){ /* part 0 pinned to A -> unordered bipartitions */ | |
| 167 | int ka=0,kb=0; | |
| 168 | for(int i=0;i<p;i++){ | |
| 169 | int inA=(i==0)?1:((m>>(i-1))&1); | |
| 170 | if(inA){ Ap[ka]=i; Af[ka]=p+i; ka++; } else { Bp[kb]=i; Bf[kb]=p+i; kb++; } | |
| 171 | } | |
| 172 | if(!ka||!kb) continue; | |
| 173 | int oa,ob; | |
| 174 | double c1=phir_mi(S,d,Ap,ka,Bf,kb,&oa); | |
| 175 | double c2=phir_mi(S,d,Bp,kb,Af,ka,&ob); | |
| 176 | if(!oa||!ob) continue; | |
| 177 | int mn=(ka<kb)?ka:kb; | |
| 178 | double cut=(c1+c2)/(double)mn; | |
| 179 | if(cut<bestcut){ bestcut=cut; bka=ka; bkb=kb; | |
| 180 | memcpy(bestA,Ap,(size_t)ka*sizeof(int)); memcpy(bestB,Bp,(size_t)kb*sizeof(int)); found=1; } | |
| 181 | } | |
| 182 | if(!found){ free(S); *ok=0; return 0.0; } | |
| 183 | int o1,o2; | |
| 184 | double I_A=phir_mi(S,d,bestA,bka,fut,p,&o1); | |
| 185 | double I_B=phir_mi(S,d,bestB,bkb,fut,p,&o2); | |
| 186 | free(S); | |
| 187 | *ok=o1&&o2; | |
| 188 | return I_whole-I_A-I_B; | |
| 189 | } | |
| 190 | ||
| 191 | /* Permutation null: shuffle the future block across rows to destroy past->future pairing. */ | |
| 192 | static void phi_r_null(const double *pairs,int n,int p,int R,double *nmean,double *nsd) | |
| 193 | { | |
| 194 | double *buf=malloc((size_t)n*2*p*sizeof(double)); | |
| 195 | double *vals=malloc((size_t)R*sizeof(double)); | |
| 196 | int *perm=malloc((size_t)n*sizeof(int)); | |
| 197 | int nv=0; | |
| 198 | for(int r=0;r<R;r++){ | |
| 199 | memcpy(buf,pairs,(size_t)n*2*p*sizeof(double)); | |
| 200 | for(int i=0;i<n;i++) perm[i]=i; | |
| 201 | for(int i=n-1;i>0;i--){ int j=(int)(phir_xs()%(unsigned long long)(i+1)); int t=perm[i]; perm[i]=perm[j]; perm[j]=t; } | |
| 202 | for(int i=0;i<n;i++) for(int k=0;k<p;k++) buf[(size_t)i*2*p+p+k]=pairs[(size_t)perm[i]*2*p+p+k]; | |
| 203 | int ok; double v=phi_r(buf,n,p,&ok); | |
| 204 | if(ok) vals[nv++]=v; | |
| 205 | } | |
| 206 | *nmean = nv? meanv(vals,nv):NAN; | |
| 207 | *nsd = nv>1? sdv(vals,nv):NAN; | |
| 208 | free(buf);free(vals);free(perm); | |
| 209 | } | |
| 210 | ||
| 211 | #endif |
| … | … | @@ -0,0 +1,542 @@ |
| 1 | /* ribophi.c — Phi_r (causal emergence) probe on ribozyme.c's replicase bootstrap. | |
| 2 | * | |
| 3 | * WHY. Pigozzi & Levin (arXiv 2607.28250) report that in the GARD model, spikes in the | |
| 4 | * PhiID causal-emergence atom Phi_r PRECEDE and predict the appearance of self-replicators, | |
| 5 | * and that Phi_r correlates with no other network/dynamical metric they tested. Our | |
| 6 | * ribozyme.c has a SHARP known ground truth: the replicase bootstrap IGNITES for compartment | |
| 7 | * size 4 and 8 (phi 6 -> ~11.5) and FAILS well-mixed and at size 16 and 32 (phi stays ~5.5-5.9). | |
| 8 | * That makes it a natural test bed: if Phi_r is a real leading indicator of an impending | |
| 9 | * transition, it should separate those arms EARLY, before the outcome is visible. | |
| 10 | * | |
| 11 | * WHAT IS COPIED VERBATIM from ribozyme.c: the RNG, phi(), replicate(), and the whole run() | |
| 12 | * dynamics (within-compartment Wright-Fisher, fidelity-from-mean-phi coupling, parasite | |
| 13 | * advantage, between-compartment selection + random assortment). Observation adds NO RNG draws, | |
| 14 | * so trajectories are bit-identical to ribozyme.c. Verified by reproducing its 3-seed table. | |
| 15 | * | |
| 16 | * MEASURE (following the paper's Materials and Methods): | |
| 17 | * Phi_r = I(X_t ; X_t+1) - sum_i I(X_t^i ; X_t+1) | |
| 18 | * i.e. Rosas/Mediano-style causal emergence: predictive information held by the WHOLE about its | |
| 19 | * own future, beyond what the parts hold about that future. The paper notes the full multivariate | |
| 20 | * estimate is intractable and uses the MINIMUM-INFORMATION BIPARTITION to reduce to two | |
| 21 | * components; we do the same (exhaustive over bipartitions here, since our part count is small): | |
| 22 | * Phi_r = I(X_t ; X_t+1) - I(A_t ; X_t+1) - I(B_t ; X_t+1) | |
| 23 | * with (A,B) the bipartition minimising the normalised cross-cut lagged MI | |
| 24 | * [ I(A_t;B_t+1) + I(B_t;A_t+1) ] / min(|A|,|B|). | |
| 25 | * Substrate: the pool-wide histogram over hairpin score phi (our analogue of GARD's "molecular | |
| 26 | * type counts" -- the compositional state), centred-log-ratio transformed with a 0.5 pseudocount, | |
| 27 | * last component dropped for full rank -- exactly the paper's preprocessing. MI is estimated | |
| 28 | * Gaussian (log-determinant form) with Ledoit-Wolf shrinkage toward the diagonal, so the | |
| 29 | * regularisation is data-driven and not a tuned knob. | |
| 30 | * | |
| 31 | * PRE-REGISTERED PREDICTIONS (written before the first run; report whatever happens): | |
| 32 | * P1 EARLY SEPARATION. Computed on the FIRST 25% of generations only (0-49 of 200), Phi_r is | |
| 33 | * HIGHER in the arms that later ignite (csz 4, 8) than in the arms that fail (well-mixed, | |
| 34 | * csz 16, 32). | |
| 35 | * P2 SPIKES PRECEDE. Igniting runs show punctuated Phi_r excursions (>3 sd above that run's own | |
| 36 | * mean) occurring BEFORE the ignition generation (first gen at which pool mean phi exceeds | |
| 37 | * its start value + 1.0). | |
| 38 | * P3 CONTROL -- the prediction that decides whether Phi_r adds anything. Raw mean phi over the | |
| 39 | * SAME first 25% does NOT separate the arms. If raw phi already separates them, Phi_r is | |
| 40 | * redundant and P1 is worthless even if it passes. This is the analogue of the paper's | |
| 41 | * "beat the raw-composition baseline" test, and it is the one that can kill the result. | |
| 42 | * | |
| 43 | * Any of these may fail. A null is a result. Do not tune to force a positive. | |
| 44 | * | |
| 45 | * build: cc -O2 -Wall -std=c11 -o ribophi ribophi.c -lm | |
| 46 | */ | |
| 47 | #include <math.h> | |
| 48 | #include <stdio.h> | |
| 49 | #include <stdint.h> | |
| 50 | #include <string.h> | |
| 51 | #include <stdlib.h> | |
| 52 | ||
| 53 | /* ---------- model constants: verbatim from ribozyme.c ---------- */ | |
| 54 | #define L 24 | |
| 55 | #define PHIMAX (L/2) /* max hairpin pairs = 12 */ | |
| 56 | #define TOTAL 4800 | |
| 57 | #define GEN 200 | |
| 58 | #define F0 0.85 | |
| 59 | #define F1 0.999 | |
| 60 | #define BETA 2.0 | |
| 61 | ||
| 62 | /* ---------- probe constants ---------- */ | |
| 63 | #define NBFULL 13 /* phi in 0..12 -> full histogram */ | |
| 64 | #define NB 8 /* merged bins: <=3,4,5,6,7,8,9,>=10 (primary substrate) */ | |
| 65 | #define P (NB-1) /* parts after dropping one CLR component */ | |
| 66 | #define PFULL (NBFULL-1) | |
| 67 | #define WIN 50 /* window length in generations = first 25% of GEN */ | |
| 68 | #define STRIDE 2 | |
| 69 | #define NSEED 20 | |
| 70 | ||
| 71 | static unsigned long long RNG; | |
| 72 | static unsigned long long xs(void){ RNG^=RNG<<13; RNG^=RNG>>7; RNG^=RNG<<17; return RNG; } | |
| 73 | static double urand(void){ return (double)(xs()>>11)*(1.0/9007199254740992.0); } | |
| 74 | ||
| 75 | 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; } | |
| 76 | 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; } | |
| 77 | ||
| 78 | static uint32_t pool[TOTAL], nxt[TOTAL]; | |
| 79 | ||
| 80 | /* ---------- the model, verbatim + observation (no RNG draws added) ---------- */ | |
| 81 | static double run(int ncomp,int csz,int group,double *f_end,double *phi_start, | |
| 82 | double histf[GEN][NBFULL], double mphi_series[GEN]) | |
| 83 | { | |
| 84 | 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; } | |
| 85 | { double sp=0; for(int i=0;i<TOTAL;i++) sp+=phi(pool[i]); *phi_start=sp/TOTAL; } | |
| 86 | ||
| 87 | for(int g=0;g<GEN;g++){ | |
| 88 | for(int cmp=0;cmp<ncomp;cmp++){ | |
| 89 | int base=cmp*csz; double sphi=0; for(int k=0;k<csz;k++) sphi+=phi(pool[base+k]); double mphi=sphi/csz; | |
| 90 | double f=F0+(F1-F0)*(mphi/PHIMAX); | |
| 91 | double tot=0; static double w[TOTAL]; static double cum[TOTAL]; | |
| 92 | 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; } | |
| 93 | 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); } | |
| 94 | } | |
| 95 | for(int i=0;i<TOTAL;i++) pool[i]=nxt[i]; | |
| 96 | ||
| 97 | if(group){ | |
| 98 | static double Pv[TOTAL/1], cumP[TOTAL/1]; double tp=0; | |
| 99 | 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; } | |
| 100 | for(int cmp=0;cmp<ncomp;cmp++){ | |
| 101 | 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; | |
| 102 | for(int k=0;k<csz;k++){ int src=par*csz+(int)(xs()%csz); nxt[cmp*csz+k]=pool[src]; } | |
| 103 | } | |
| 104 | for(int i=0;i<TOTAL;i++) pool[i]=nxt[i]; | |
| 105 | } | |
| 106 | ||
| 107 | /* --- observation only: pool-wide phi histogram + mean phi. No RNG consumed. --- */ | |
| 108 | if(histf){ | |
| 109 | for(int b=0;b<NBFULL;b++) histf[g][b]=0.0; | |
| 110 | double sp=0; | |
| 111 | for(int i=0;i<TOTAL;i++){ int p=phi(pool[i]); histf[g][p]+=1.0; sp+=p; } | |
| 112 | if(mphi_series) mphi_series[g]=sp/TOTAL; | |
| 113 | } | |
| 114 | } | |
| 115 | double sp=0; for(int i=0;i<TOTAL;i++) sp+=phi(pool[i]); double mphi=sp/TOTAL; | |
| 116 | *f_end=F0+(F1-F0)*(mphi/PHIMAX); return mphi; | |
| 117 | } | |
| 118 | ||
| 119 | /* ---------- compositional preprocessing: merge bins, CLR, drop last ---------- */ | |
| 120 | static void merge_bins(const double *full, double *out) /* NBFULL -> NB */ | |
| 121 | { | |
| 122 | for(int b=0;b<NB;b++) out[b]=0.0; | |
| 123 | for(int p=0;p<NBFULL;p++){ | |
| 124 | int b; | |
| 125 | if(p<=3) b=0; else if(p>=10) b=NB-1; else b=p-3; /* 4,5,6,7,8,9 -> 1..6 */ | |
| 126 | out[b]+=full[p]; | |
| 127 | } | |
| 128 | } | |
| 129 | static void clr(const double *cnt,int nb,double *out) /* -> nb-1 components */ | |
| 130 | { | |
| 131 | double tot=0; for(int i=0;i<nb;i++) tot+=cnt[i]+0.5; | |
| 132 | double lg[NBFULL]; double mean=0; | |
| 133 | for(int i=0;i<nb;i++){ lg[i]=log((cnt[i]+0.5)/tot); mean+=lg[i]; } | |
| 134 | mean/=nb; | |
| 135 | for(int i=0;i<nb-1;i++) out[i]=lg[i]-mean; | |
| 136 | } | |
| 137 | ||
| 138 | /* ---------- Gaussian MI via log-determinants, Ledoit-Wolf shrinkage ---------- */ | |
| 139 | /* Z: n x d row-major (already the joint [past | future] vector). Builds shrunk covariance. */ | |
| 140 | static void covariance_lw(const double *Z,int n,int d,double *S) | |
| 141 | { | |
| 142 | double *mu=calloc(d,sizeof(double)); | |
| 143 | for(int t=0;t<n;t++) for(int j=0;j<d;j++) mu[j]+=Z[t*d+j]; | |
| 144 | for(int j=0;j<d;j++) mu[j]/=n; | |
| 145 | double *Y=malloc((size_t)n*d*sizeof(double)); | |
| 146 | for(int t=0;t<n;t++) for(int j=0;j<d;j++) Y[t*d+j]=Z[t*d+j]-mu[j]; | |
| 147 | ||
| 148 | for(int i=0;i<d;i++) for(int j=0;j<d;j++){ | |
| 149 | double s=0; for(int t=0;t<n;t++) s+=Y[t*d+i]*Y[t*d+j]; | |
| 150 | S[i*d+j]=s/n; | |
| 151 | } | |
| 152 | /* Ledoit-Wolf shrinkage intensity toward the diagonal target */ | |
| 153 | double num=0, den=0; | |
| 154 | for(int i=0;i<d;i++) for(int j=0;j<d;j++){ | |
| 155 | if(i==j) continue; | |
| 156 | double sij=S[i*d+j], v=0; | |
| 157 | for(int t=0;t<n;t++){ double e=Y[t*d+i]*Y[t*d+j]-sij; v+=e*e; } | |
| 158 | v*= (double)n/((double)(n-1)*(n-1)*(n-1)); | |
| 159 | num+=v; den+=sij*sij; | |
| 160 | } | |
| 161 | double lam = (den>0)? num/den : 1.0; | |
| 162 | if(lam<0) lam=0; | |
| 163 | if(lam>1) lam=1; | |
| 164 | for(int i=0;i<d;i++) for(int j=0;j<d;j++) if(i!=j) S[i*d+j]*=(1.0-lam); | |
| 165 | /* unbiased-ish rescale to sample covariance convention */ | |
| 166 | double sc=(double)n/(double)(n-1); | |
| 167 | for(int i=0;i<d;i++) for(int j=0;j<d;j++) S[i*d+j]*=sc; | |
| 168 | free(mu); free(Y); | |
| 169 | } | |
| 170 | /* log|S[idx,idx]| via Cholesky, with jitter fallback. Returns 0 and sets *ok=0 on failure. */ | |
| 171 | static double logdet_sub(const double *S,int d,const int *idx,int k,int *ok) | |
| 172 | { | |
| 173 | double *M=malloc((size_t)k*k*sizeof(double)); | |
| 174 | for(int a=0;a<k;a++) for(int b=0;b<k;b++) M[a*k+b]=S[idx[a]*d+idx[b]]; | |
| 175 | double tr=0; for(int a=0;a<k;a++) tr+=M[a*k+a]; | |
| 176 | double jitter=1e-10*(tr/k+1e-12); | |
| 177 | for(int attempt=0;attempt<8;attempt++){ | |
| 178 | double *C=malloc((size_t)k*k*sizeof(double)); | |
| 179 | memcpy(C,M,(size_t)k*k*sizeof(double)); | |
| 180 | for(int a=0;a<k;a++) C[a*k+a]+=jitter; | |
| 181 | int bad=0; double ld=0; | |
| 182 | for(int a=0;a<k && !bad;a++){ | |
| 183 | for(int b=0;b<=a;b++){ | |
| 184 | double s=C[a*k+b]; | |
| 185 | for(int m=0;m<b;m++) s-=C[a*k+m]*C[b*k+m]; | |
| 186 | if(a==b){ if(s<=0){ bad=1; break; } C[a*k+a]=sqrt(s); ld+=2.0*log(C[a*k+a]); } | |
| 187 | else C[a*k+b]=s/C[b*k+b]; | |
| 188 | } | |
| 189 | } | |
| 190 | free(C); | |
| 191 | if(!bad){ free(M); *ok=1; return ld; } | |
| 192 | jitter = (jitter>0)? jitter*100.0 : 1e-9*(tr/k+1e-12); | |
| 193 | } | |
| 194 | free(M); *ok=0; return 0.0; | |
| 195 | } | |
| 196 | static double mi_sub(const double *S,int d,const int *A,int ka,const int *B,int kb,int *ok) | |
| 197 | { | |
| 198 | int *U=malloc((size_t)(ka+kb)*sizeof(int)); | |
| 199 | memcpy(U,A,(size_t)ka*sizeof(int)); memcpy(U+ka,B,(size_t)kb*sizeof(int)); | |
| 200 | int o1,o2,o3; | |
| 201 | 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); | |
| 202 | free(U); | |
| 203 | *ok = o1&&o2&&o3; | |
| 204 | return 0.5*(la+lb-lu); | |
| 205 | } | |
| 206 | ||
| 207 | /* Phi_r over n transition pairs of p-dim CLR vectors. | |
| 208 | pairs: n x (2p) row-major, [past(0..p-1) | future(p..2p-1)]. Returns 0 and *ok=0 on failure. */ | |
| 209 | static double phi_r(const double *pairs,int n,int p,int *ok) | |
| 210 | { | |
| 211 | int d=2*p; | |
| 212 | double *S=malloc((size_t)d*d*sizeof(double)); | |
| 213 | covariance_lw(pairs,n,d,S); | |
| 214 | ||
| 215 | int *past=malloc(p*sizeof(int)), *fut=malloc(p*sizeof(int)); | |
| 216 | for(int i=0;i<p;i++){ past[i]=i; fut[i]=p+i; } | |
| 217 | ||
| 218 | int o; | |
| 219 | double I_whole = mi_sub(S,d,past,p,fut,p,&o); | |
| 220 | if(!o){ free(S);free(past);free(fut); *ok=0; return 0.0; } | |
| 221 | ||
| 222 | /* minimum-information bipartition over the p parts */ | |
| 223 | int bestA[P+1], bestB[P+1], bka=0, bkb=0; double bestcut=1e300; int found=0; | |
| 224 | int *Ap=malloc(p*sizeof(int)), *Bp=malloc(p*sizeof(int)); | |
| 225 | int *Af=malloc(p*sizeof(int)), *Bf=malloc(p*sizeof(int)); | |
| 226 | for(unsigned m=1; m < (1u<<(p-1)); m++){ /* fix part 0 in A -> unordered bipartitions */ | |
| 227 | int ka=0,kb=0; | |
| 228 | for(int i=0;i<p;i++){ | |
| 229 | int inA = (i==0)? 1 : ((m>>(i-1))&1); | |
| 230 | if(inA){ Ap[ka]=i; Af[ka]=p+i; ka++; } else { Bp[kb]=i; Bf[kb]=p+i; kb++; } | |
| 231 | } | |
| 232 | if(ka==0||kb==0) continue; | |
| 233 | int oa,ob; | |
| 234 | double c1=mi_sub(S,d,Ap,ka,Bf,kb,&oa); | |
| 235 | double c2=mi_sub(S,d,Bp,kb,Af,ka,&ob); | |
| 236 | if(!oa||!ob) continue; | |
| 237 | int mn = (ka<kb)?ka:kb; | |
| 238 | double cut=(c1+c2)/(double)mn; | |
| 239 | if(cut<bestcut){ bestcut=cut; bka=ka; bkb=kb; | |
| 240 | memcpy(bestA,Ap,(size_t)ka*sizeof(int)); memcpy(bestB,Bp,(size_t)kb*sizeof(int)); found=1; } | |
| 241 | } | |
| 242 | free(Ap);free(Bp);free(Af);free(Bf); | |
| 243 | if(!found){ free(S);free(past);free(fut); *ok=0; return 0.0; } | |
| 244 | ||
| 245 | int o1,o2; | |
| 246 | double I_A = mi_sub(S,d,bestA,bka,fut,p,&o1); | |
| 247 | double I_B = mi_sub(S,d,bestB,bkb,fut,p,&o2); | |
| 248 | free(S);free(past);free(fut); | |
| 249 | *ok = o1&&o2; | |
| 250 | return I_whole - I_A - I_B; | |
| 251 | } | |
| 252 | ||
| 253 | /* ---------- helpers ---------- */ | |
| 254 | static int cmpd(const void*a,const void*b){ double x=*(const double*)a,y=*(const double*)b; return (x>y)-(x<y); } | |
| 255 | 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; } | |
| 256 | static double meanv(const double *v,int n){ double s=0; for(int i=0;i<n;i++) s+=v[i]; return s/n; } | |
| 257 | 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)); } | |
| 258 | /* Mann-Whitney U -> normal approximation two-sided p (no ties correction; n=20 per group) */ | |
| 259 | static double mannwhitney_p(const double *x,int nx,const double *y,int ny) | |
| 260 | { | |
| 261 | double U=0; | |
| 262 | 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; } | |
| 263 | double mu=0.5*nx*ny, sd=sqrt((double)nx*ny*(nx+ny+1)/12.0); | |
| 264 | if(sd<=0) return 1.0; | |
| 265 | double z=fabs(U-mu)/sd; | |
| 266 | return erfc(z/sqrt(2.0)); | |
| 267 | } | |
| 268 | ||
| 269 | typedef struct { const char*name; int ncomp,csz,group; } Arm; | |
| 270 | ||
| 271 | /* Permutation null: destroy the past->future pairing by shuffling the future block across rows, | |
| 272 | recompute Phi_r R times, return mean and sd of the null distribution. */ | |
| 273 | static void phi_r_null(const double *pairs,int n,int p,int R,double *nmean,double *nsd) | |
| 274 | { | |
| 275 | double *buf=malloc((size_t)n*2*p*sizeof(double)); | |
| 276 | double *vals=malloc((size_t)R*sizeof(double)); | |
| 277 | int *perm=malloc((size_t)n*sizeof(int)); | |
| 278 | int nv=0; | |
| 279 | for(int r=0;r<R;r++){ | |
| 280 | memcpy(buf,pairs,(size_t)n*2*p*sizeof(double)); | |
| 281 | for(int i=0;i<n;i++) perm[i]=i; | |
| 282 | 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; } | |
| 283 | 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]; | |
| 284 | int ok; double v=phi_r(buf,n,p,&ok); | |
| 285 | if(ok) vals[nv++]=v; | |
| 286 | } | |
| 287 | *nmean = nv? meanv(vals,nv):NAN; | |
| 288 | *nsd = nv>1? sdv(vals,nv):NAN; | |
| 289 | free(buf);free(vals);free(perm); | |
| 290 | } | |
| 291 | ||
| 292 | int main(void) | |
| 293 | { | |
| 294 | static Arm arms[]={ | |
| 295 | {"well-mixed", 1, TOTAL, 0}, | |
| 296 | {"comp 1200x4", 1200, 4, 1}, | |
| 297 | {"comp 600x8", 600, 8, 1}, | |
| 298 | {"comp 300x16", 300, 16, 1}, | |
| 299 | {"comp 150x32", 150, 32, 1}, | |
| 300 | }; | |
| 301 | const int NA=(int)(sizeof(arms)/sizeof(arms[0])); | |
| 302 | ||
| 303 | static double histf[GEN][NBFULL]; | |
| 304 | static double mser[GEN]; | |
| 305 | ||
| 306 | /* ============ 0a. estimator validation on synthetic data with KNOWN answers ============ | |
| 307 | A null result downstream is uninterpretable unless the estimator is shown to respond | |
| 308 | with the right SIGN and SENSITIVITY first. Three systems, same p and n as the real probe: | |
| 309 | (a) independent -> Phi_r ~ 0 | |
| 310 | (b) common driver -> Phi_r < 0 (redundancy dominates: each half already predicts) | |
| 311 | (c) difference -> Phi_r > 0 (synergy: neither half alone reveals d1-d2) */ | |
| 312 | printf("=== 0a. VALIDATE estimator on synthetic systems with known answers (p=%d, n=100) ===\n",P); | |
| 313 | { | |
| 314 | RNG=0xabcdef0123456789ULL; | |
| 315 | int n=100; static double sy[100][2*P]; | |
| 316 | double gz[4]; | |
| 317 | /* (a) independent */ | |
| 318 | 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); } | |
| 319 | int ok; double va=phi_r(&sy[0][0],n,P,&ok); | |
| 320 | /* (b) common driver: every part = z + small noise, at both times */ | |
| 321 | for(int t=0;t<n;t++){ | |
| 322 | double u=urand(),v=urand(); if(u<1e-300)u=1e-300; double z=sqrt(-2*log(u))*cos(6.283185307*v); | |
| 323 | 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; } | |
| 324 | } | |
| 325 | int okb; double vb=phi_r(&sy[0][0],n,P,&okb); | |
| 326 | /* (c) synergy under EVERY bipartition. Each part is z + 0.1*d_i, so z swamps it. | |
| 327 | future_0 = sum_i w_i d_i with weights w summing to ZERO, so the contrast is invariant | |
| 328 | to the common offset z: the WHOLE recovers it exactly (z cancels), while any proper | |
| 329 | subset has sum(w) != 0 in general and cannot. Unlike a single synergistic pair, this | |
| 330 | cannot be defeated by the MIB choosing a cut that keeps the pair together. */ | |
| 331 | { | |
| 332 | static const double w[P]={+1,-1,+1,-1,+1,-1,0}; | |
| 333 | double dd[P]; | |
| 334 | for(int t=0;t<n;t++){ | |
| 335 | double u=urand(),v=urand(); if(u<1e-300)u=1e-300; double z=sqrt(-2*log(u))*cos(6.283185307*v); | |
| 336 | 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); } | |
| 337 | for(int i=0;i<P;i++) sy[t][i] = z + 0.1*dd[i]; | |
| 338 | 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; } | |
| 339 | double acc=0; for(int i=0;i<P;i++) acc+=w[i]*dd[i]; | |
| 340 | sy[t][P+0]=acc; | |
| 341 | } | |
| 342 | } | |
| 343 | int okc; double vc=phi_r(&sy[0][0],n,P,&okc); | |
| 344 | (void)gz; | |
| 345 | /* judge each against its OWN permutation null rather than an arbitrary absolute cut */ | |
| 346 | double nma,nsa,nmb,nsb,nmc,nsc; | |
| 347 | phi_r_null(&sy[0][0],n,P,40,&nmc,&nsc); /* sy currently holds (c) */ | |
| 348 | double zc = (nsc>0)? (vc-nmc)/nsc : NAN; | |
| 349 | printf(" (c) synergy (zero-sum contrast) -> Phi_r = %+8.5f vs shuffled null %+.5f +- %.5f z=%+.1f\n", vc,nmc,nsc,zc); | |
| 350 | /* regenerate (a) and (b) to get their nulls too */ | |
| 351 | 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); } | |
| 352 | phi_r_null(&sy[0][0],n,P,40,&nma,&nsa); | |
| 353 | double za=(nsa>0)?(va-nma)/nsa:NAN; | |
| 354 | for(int t=0;t<n;t++){ | |
| 355 | double u=urand(),v2=urand(); if(u<1e-300)u=1e-300; double z=sqrt(-2*log(u))*cos(6.283185307*v2); | |
| 356 | 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; } | |
| 357 | } | |
| 358 | phi_r_null(&sy[0][0],n,P,40,&nmb,&nsb); | |
| 359 | double zb=(nsb>0)?(vb-nmb)/nsb:NAN; | |
| 360 | printf(" (a) independent parts -> Phi_r = %+8.5f vs shuffled null %+.5f +- %.5f z=%+.1f expect ~0\n", va,nma,nsa,za); | |
| 361 | printf(" (b) common driver -> Phi_r = %+8.5f vs shuffled null %+.5f +- %.5f z=%+.1f expect <0\n", vb,nmb,nsb,zb); | |
| 362 | (void)ok;(void)okb;(void)okc; | |
| 363 | int pa=(fabs(za)<3.0), pb=(zb<-3.0), pc=(zc>3.0); | |
| 364 | printf(" verdict: (a) %s (b) %s (c) %s\n", pa?"PASS":"FAIL", pb?"PASS":"FAIL", pc?"PASS":"FAIL"); | |
| 365 | if(pa&&pb&&pc) printf(" => estimator recovers all three signs against its own null. A null below is informative.\n"); | |
| 366 | else printf(" => ESTIMATOR FAILS A SIGN TEST. Nulls below are confounded and must NOT be read as\n claims about Phi_r itself.\n"); | |
| 367 | 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"); | |
| 368 | } | |
| 369 | ||
| 370 | /* ================= 0. reproduce ribozyme.c on its own 3 seeds ================= */ | |
| 371 | printf("=== 0. VERIFY: instrumented model reproduces ribozyme.c (same 3 seeds) ===\n"); | |
| 372 | printf(" arm | phi start -> end | fidelity end (ribozyme.c reference)\n"); | |
| 373 | unsigned long long vs[]={0x9e3779b97f4a7c15ULL,0xd1b54a32d192ed03ULL,0x2545f4914f6cdd1dULL}; | |
| 374 | for(int a=0;a<NA;a++){ | |
| 375 | double S=0,E=0,F=0; | |
| 376 | for(int s=0;s<3;s++){ double fe,ps; RNG=vs[s]; | |
| 377 | E+=run(arms[a].ncomp,arms[a].csz,arms[a].group,&fe,&ps,histf,mser); S+=ps; F+=fe; } | |
| 378 | printf(" %-15s | %.2f -> %5.2f | %.3f\n",arms[a].name,S/3,E/3,F/3); | |
| 379 | } | |
| 380 | ||
| 381 | /* ================= 1. run NSEED seeds per arm, collect Phi_r ================= */ | |
| 382 | printf("\n=== 1. Phi_r probe: %d seeds x %d arms, GEN=%d, window=%d gens (first 25%%) ===\n",NSEED,NA,GEN,WIN); | |
| 383 | printf(" substrate = pool-wide phi histogram (%d merged bins -> %d CLR parts), Gaussian MI,\n",NB,P); | |
| 384 | printf(" Ledoit-Wolf shrinkage, exhaustive minimum-information bipartition.\n\n"); | |
| 385 | ||
| 386 | static double earlyPhi[8][NSEED]; /* Phi_r on gens 0..WIN-1 */ | |
| 387 | static double earlyRawPhi[8][NSEED]; /* CONTROL: raw mean phi, same gens */ | |
| 388 | static double finalPhi[8][NSEED]; /* outcome: end mean phi */ | |
| 389 | static int spikeBefore[8], spikeAny[8], ignited[8]; | |
| 390 | ||
| 391 | static int ignGen[8][NSEED]; /* generation at which each run ignites (-1 = never) */ | |
| 392 | static double clrAll[8][NSEED][GEN][P];/* CLR series kept for the pre-ignition test */ | |
| 393 | static double mserAll[8][NSEED][GEN]; | |
| 394 | ||
| 395 | for(int a=0;a<NA;a++){ | |
| 396 | spikeBefore[a]=spikeAny[a]=ignited[a]=0; | |
| 397 | for(int s=0;s<NSEED;s++){ | |
| 398 | double fe,ps; | |
| 399 | RNG = 0x9e3779b97f4a7c15ULL ^ (0x100000001b3ULL*(unsigned long long)(s+1)); | |
| 400 | for(int k=0;k<3;k++) xs(); /* decorrelate seed streams */ | |
| 401 | double endphi = run(arms[a].ncomp,arms[a].csz,arms[a].group,&fe,&ps,histf,mser); | |
| 402 | finalPhi[a][s]=endphi; | |
| 403 | ||
| 404 | /* CONTROL: raw compositional signal over the same early window */ | |
| 405 | double rp=0; for(int g=0;g<WIN;g++) rp+=mser[g]; | |
| 406 | earlyRawPhi[a][s]=rp/WIN; | |
| 407 | ||
| 408 | /* build CLR series once */ | |
| 409 | static double clrser[GEN][P]; | |
| 410 | for(int g=0;g<GEN;g++){ double m8[NB]; merge_bins(histf[g],m8); clr(m8,NB,clrser[g]); } | |
| 411 | memcpy(clrAll[a][s],clrser,sizeof(clrser)); | |
| 412 | memcpy(mserAll[a][s],mser,sizeof(mser)); | |
| 413 | { int gg=-1; for(int g=0;g<GEN;g++) if(mser[g] > ps+1.0){ gg=g; break; } ignGen[a][s]=gg; } | |
| 414 | ||
| 415 | /* Phi_r on the early window only */ | |
| 416 | { | |
| 417 | int n=WIN-1; static double pairs[GEN][2*P]; | |
| 418 | 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]; } } | |
| 419 | int ok; double v=phi_r(&pairs[0][0],n,P,&ok); | |
| 420 | earlyPhi[a][s] = ok? v : NAN; | |
| 421 | } | |
| 422 | ||
| 423 | /* full sliding trajectory, for the spike test */ | |
| 424 | { | |
| 425 | int nw=0; static double traj[GEN]; static int tcenter[GEN]; | |
| 426 | for(int st=0; st+WIN<=GEN; st+=STRIDE){ | |
| 427 | int n=WIN-1; static double pairs[GEN][2*P]; | |
| 428 | 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]; } } | |
| 429 | int ok; double v=phi_r(&pairs[0][0],n,P,&ok); | |
| 430 | if(ok){ traj[nw]=v; tcenter[nw]=st+WIN/2; nw++; } | |
| 431 | } | |
| 432 | (void)meanv; | |
| 433 | /* ignition generation: first gen where mean phi exceeds start + 1.0 */ | |
| 434 | int gign=-1; for(int g=0;g<GEN;g++) if(mser[g] > ps+1.0){ gign=g; break; } | |
| 435 | if(gign>=0) ignited[a]++; | |
| 436 | if(nw>3){ | |
| 437 | double mu=meanv(traj,nw), sd=sdv(traj,nw); | |
| 438 | for(int i=0;i<nw;i++) if(traj[i] > mu+3.0*sd){ | |
| 439 | spikeAny[a]++; | |
| 440 | if(gign>=0 && tcenter[i]<gign) spikeBefore[a]++; | |
| 441 | break; | |
| 442 | } | |
| 443 | } | |
| 444 | } | |
| 445 | } | |
| 446 | } | |
| 447 | ||
| 448 | /* ================= 2. results ================= */ | |
| 449 | printf(" arm | end phi | P3 CONTROL raw phi | P1 Phi_r (early, gens 0-%d)\n",WIN-1); | |
| 450 | printf(" | (gen%3d)| mean+-sd (0-%2d) | median mean+-sd\n",GEN,WIN-1); | |
| 451 | printf(" ----------------+---------+--------------------+---------------------------\n"); | |
| 452 | for(int a=0;a<NA;a++){ | |
| 453 | double ep[NSEED]; int ne=0; | |
| 454 | for(int s=0;s<NSEED;s++) if(!isnan(earlyPhi[a][s])) ep[ne++]=earlyPhi[a][s]; | |
| 455 | printf(" %-15s | %6.2f | %5.2f +- %4.2f | %7.4f %7.4f +- %6.4f (n=%d)\n", | |
| 456 | arms[a].name, meanv(finalPhi[a],NSEED), | |
| 457 | meanv(earlyRawPhi[a],NSEED), sdv(earlyRawPhi[a],NSEED), | |
| 458 | ne?median(ep,ne):NAN, ne?meanv(ep,ne):NAN, ne?sdv(ep,ne):NAN, ne); | |
| 459 | } | |
| 460 | ||
| 461 | /* pool the igniting arms vs the failing arms */ | |
| 462 | double ig[2*NSEED], fl[3*NSEED]; int nig=0,nfl=0; | |
| 463 | double igR[2*NSEED], flR[3*NSEED]; int nigR=0,nflR=0; | |
| 464 | for(int a=0;a<NA;a++){ | |
| 465 | int isIgnite = (arms[a].csz==4 || arms[a].csz==8) && arms[a].group; | |
| 466 | for(int s=0;s<NSEED;s++){ | |
| 467 | if(!isnan(earlyPhi[a][s])){ if(isIgnite) ig[nig++]=earlyPhi[a][s]; else fl[nfl++]=earlyPhi[a][s]; } | |
| 468 | if(isIgnite) igR[nigR++]=earlyRawPhi[a][s]; else flR[nflR++]=earlyRawPhi[a][s]; | |
| 469 | } | |
| 470 | } | |
| 471 | printf("\n --- P1: early Phi_r, igniting arms (csz 4,8) vs failing arms (well-mixed,16,32) ---\n"); | |
| 472 | printf(" igniting: mean %.4f (n=%d) failing: mean %.4f (n=%d) Mann-Whitney p=%.4g\n", | |
| 473 | meanv(ig,nig),nig, meanv(fl,nfl),nfl, mannwhitney_p(ig,nig,fl,nfl)); | |
| 474 | printf(" --- P3 CONTROL: raw mean phi over the SAME early window ---\n"); | |
| 475 | printf(" igniting: mean %.4f (n=%d) failing: mean %.4f (n=%d) Mann-Whitney p=%.4g\n", | |
| 476 | meanv(igR,nigR),nigR, meanv(flR,nflR),nflR, mannwhitney_p(igR,nigR,flR,nflR)); | |
| 477 | ||
| 478 | printf("\n --- P2: spikes (>3 sd above run mean) and whether they precede ignition ---\n"); | |
| 479 | printf(" arm | runs ignited | runs with a spike | spike BEFORE ignition\n"); | |
| 480 | for(int a=0;a<NA;a++) | |
| 481 | printf(" %-15s | %5d / %-5d | %6d / %-8d | %d\n", | |
| 482 | arms[a].name, ignited[a],NSEED, spikeAny[a],NSEED, spikeBefore[a]); | |
| 483 | ||
| 484 | /* ================= 3. WHEN does ignition happen? ================= */ | |
| 485 | printf("\n=== 3. Ignition timing (first gen where pool mean phi > start+1.0) ===\n"); | |
| 486 | printf(" arm | ignited | ignition gen: min / median / max\n"); | |
| 487 | int gmin_global=GEN; | |
| 488 | for(int a=0;a<NA;a++){ | |
| 489 | double gv[NSEED]; int ng=0; int mn=GEN,mx=-1; | |
| 490 | 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]; } | |
| 491 | if(ng){ printf(" %-15s | %2d/%2d | %d / %.1f / %d\n",arms[a].name,ng,NSEED,mn,median(gv,ng),mx); | |
| 492 | if(mn<gmin_global) gmin_global=mn; } | |
| 493 | else printf(" %-15s | %2d/%2d | -\n",arms[a].name,ng,NSEED); | |
| 494 | } | |
| 495 | printf(" => earliest ignition anywhere: gen %d. The window used in section 1 was gens 0-%d,\n",gmin_global,WIN-1); | |
| 496 | printf(" so that window was NOT pre-transition. Re-testing strictly before gen %d.\n",gmin_global); | |
| 497 | ||
| 498 | /* ================= 4. strictly PRE-IGNITION window ================= */ | |
| 499 | int PRE = gmin_global; /* gens 0..PRE-1 contain no ignited state in ANY run */ | |
| 500 | printf("\n=== 4. P1 re-test on a strictly PRE-IGNITION window (gens 0-%d) ===\n",PRE-1); | |
| 501 | if(PRE < 4){ | |
| 502 | printf(" Only %d pre-ignition generations exist -- too few transitions per run for a per-run\n",PRE); | |
| 503 | printf(" estimate. Pooling transition pairs across the %d seeds of each arm instead.\n",NSEED); | |
| 504 | } | |
| 505 | { | |
| 506 | int perRun = PRE-1; /* transition pairs available per run */ | |
| 507 | if(perRun < 1){ | |
| 508 | printf(" NO pre-ignition transitions exist at generation resolution. Test cannot be run here.\n"); | |
| 509 | } else { | |
| 510 | printf(" arm | pooled n | Phi_r (pre-ignition) | CONTROL raw phi mean+-sd\n"); | |
| 511 | printf(" ----------------+----------+----------------------+-------------------------\n"); | |
| 512 | static double poolPairs[NSEED*GEN][2*P]; | |
| 513 | double preRawIg[NSEED*8]; int nRawIg=0; double preRawFl[NSEED*8]; int nRawFl=0; | |
| 514 | double phirIg[8]; double phirFl[8]; int nIgA=0,nFlA=0; | |
| 515 | for(int a=0;a<NA;a++){ | |
| 516 | int n=0; | |
| 517 | for(int s=0;s<NSEED;s++) | |
| 518 | for(int t=0;t<perRun;t++){ | |
| 519 | 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]; } | |
| 520 | n++; | |
| 521 | } | |
| 522 | int ok; double v=phi_r(&poolPairs[0][0],n,P,&ok); | |
| 523 | double nm,ns; phi_r_null(&poolPairs[0][0],n,P,40,&nm,&ns); | |
| 524 | double zz=(ns>0)?(v-nm)/ns:NAN; | |
| 525 | double rr[NSEED]; | |
| 526 | 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; } | |
| 527 | printf(" %-15s | %8d | %9.5f (z=%+5.1f) | %.3f +- %.3f\n", | |
| 528 | arms[a].name,n, ok? v : NAN, zz, meanv(rr,NSEED), sdv(rr,NSEED)); | |
| 529 | int isIgnite = (arms[a].csz==4 || arms[a].csz==8) && arms[a].group; | |
| 530 | if(isIgnite){ if(ok) phirIg[nIgA++]=v; for(int s=0;s<NSEED;s++) preRawIg[nRawIg++]=rr[s]; } | |
| 531 | else { if(ok) phirFl[nFlA++]=v; for(int s=0;s<NSEED;s++) preRawFl[nRawFl++]=rr[s]; } | |
| 532 | } | |
| 533 | printf("\n P1 (pre-ignition): mean Phi_r igniting arms %.5f vs failing arms %.5f\n", | |
| 534 | meanv(phirIg,nIgA), meanv(phirFl,nFlA)); | |
| 535 | printf(" P3 CONTROL (pre-ignition): raw phi igniting %.4f vs failing %.4f Mann-Whitney p=%.4g\n", | |
| 536 | meanv(preRawIg,nRawIg), meanv(preRawFl,nRawFl), mannwhitney_p(preRawIg,nRawIg,preRawFl,nRawFl)); | |
| 537 | } | |
| 538 | } | |
| 539 | ||
| 540 | printf("\n=== verdict is read off the pre-registered predictions above; see ribophi.md ===\n"); | |
| 541 | return 0; | |
| 542 | } |
| … | … | @@ -0,0 +1,88 @@ |
| 1 | # Φ_r probe against `ribozyme.c` — honest null, with a diagnosis | |
| 2 | ||
| 3 | **Code:** `ribophi.c` · **Prompted by:** Pigozzi & Levin, *Causal Architecture Dynamics Prior to | |
| 4 | Arrival of Self-replicators in a Model of Catalytic Networks Relevant to Origin-of-Life* | |
| 5 | (arXiv 2607.28250, q-bio.PE). | |
| 6 | ||
| 7 | ## What was tested | |
| 8 | ||
| 9 | The paper reports that in GARD, spikes in the ΦID causal-emergence atom Φ_r **precede and predict** | |
| 10 | the appearance of self-replicators, and that Φ_r beats raw-composition baselines at that prediction. | |
| 11 | `ribozyme.c` has a sharp known ground truth — the replicase bootstrap ignites at compartment size 4 | |
| 12 | and 8 and fails well-mixed / 16 / 32 — so it is a natural test bed for a claimed leading indicator. | |
| 13 | ||
| 14 | Measure, following the paper's Materials and Methods: Φ_r = I(X_t;X_t+1) − Σᵢ I(X_t^i;X_t+1), | |
| 15 | reduced via the minimum-information bipartition to Φ_r = I(X_t;X_t+1) − I(A_t;X_t+1) − I(B_t;X_t+1). | |
| 16 | Substrate = pool-wide histogram over hairpin score φ (our analogue of GARD's molecular type counts), | |
| 17 | CLR-transformed with a 0.5 pseudocount, last component dropped for full rank. Gaussian MI via | |
| 18 | log-determinants, Ledoit-Wolf shrinkage (data-driven, not a tuned knob), exhaustive bipartition. | |
| 19 | ||
| 20 | Model dynamics copied **verbatim**; observation consumes no RNG, so trajectories are bit-identical. | |
| 21 | Confirmed by reproducing `ribozyme.c`'s 3-seed table exactly: 5.50 / 11.64 / 11.34 / 5.90 / 5.73. | |
| 22 | ||
| 23 | ## Estimator validation (before reading any null) | |
| 24 | ||
| 25 | | synthetic system | Φ_r | vs shuffled null | verdict | | |
| 26 | |---|---|---|---| | |
| 27 | | independent parts | +0.00000 | z = +0.2 | PASS (≈0) | | |
| 28 | | common driver (redundancy) | −1.78450 | z = −165.3 | PASS (<0) | | |
| 29 | | zero-sum contrast (synergy) | +0.03484 | z = +3.4 | PASS (>0) | | |
| 30 | ||
| 31 | All three signs recovered. **Caveat carried forward:** |z| for redundancy exceeds that for synergy by | |
| 32 | ~50×. This Gaussian estimator sees redundancy far more readily than synergy. | |
| 33 | ||
| 34 | ## Results — all three pre-registered predictions | |
| 35 | ||
| 36 | **P1 (early separation) — FAILS.** On a strictly pre-ignition window (gens 0–5, before *any* run has | |
| 37 | ignited), Φ_r does not order by outcome: | |
| 38 | ||
| 39 | | arm | Φ_r pre-ignition (z vs shuffled) | outcome | | |
| 40 | |---|---|---| | |
| 41 | | well-mixed | −0.06979 (z=−3.6) | fails | | |
| 42 | | comp 1200×4 | −0.02660 (z=−1.5) | **ignites** | | |
| 43 | | comp 600×8 | −0.06709 (z=−4.2) | **ignites** | | |
| 44 | | comp 300×16 | −0.03239 (z=−9.6) | fails | | |
| 45 | | comp 150×32 | +0.00048 (z=+0.0) | fails | | |
| 46 | ||
| 47 | Igniting arms mean −0.047 vs failing −0.034 — the *wrong direction*, and the ordering is scrambled | |
| 48 | (csz=4 sits between two failing arms; csz=8 sits next to well-mixed). The z-scores show real | |
| 49 | past→future structure exists, but the arm with the strongest signal (csz=16, z=−9.6) is one that fails. | |
| 50 | ||
| 51 | **P2 (spikes precede) — FAILS.** 0 of 100 runs had a Φ_r spike (>3 sd above run mean) preceding | |
| 52 | ignition. Only 8/100 runs had a spike at all. | |
| 53 | ||
| 54 | **P3 (control) — FAILS, and this is the one that matters.** Raw mean φ over the same window already | |
| 55 | separates the arms at p=3.1e-17, against Φ_r's p=1.1e-6 in the first-25% window. Raw composition | |
| 56 | beats Φ_r by eleven orders of magnitude in p. Φ_r adds nothing here. | |
| 57 | ||
| 58 | ## Diagnosis — this is about the model, not about Φ_r | |
| 59 | ||
| 60 | Ignition timing, measured: **csz=4 ignites at gen 6–10 (median 7); csz=8 at gen 12–151 (median 48.5).** | |
| 61 | The paper's protocol — take the first 25% of the trajectory and predict the rest — assumes a | |
| 62 | *quiescent pre-replicative phase*. `ribozyme.c` has none. The first 25% (gens 0–49) is already | |
| 63 | post-transition for csz=4. Selection is on from generation 1, and the arms diverge immediately: even | |
| 64 | in gens 0–5 raw φ differs (6.46 vs 5.62). | |
| 65 | ||
| 66 | So `ribozyme.c` is the **wrong test bed** for a leading-indicator claim. It does not falsify the | |
| 67 | paper; it fails to host the test. | |
| 68 | ||
| 69 | ## The one positive finding, with the right sign flagged | |
| 70 | ||
| 71 | *During* the transition (gens 0–49), Φ_r does respond strongly — igniting −0.217 vs failing −0.033, | |
| 72 | p=1.1e-6. But it moves **more negative**, i.e. redundancy-dominated. In `ribozyme.c` the transition | |
| 73 | registers as the population *homogenising* onto high-φ folds, not as a rise in synergy. The paper | |
| 74 | reports Φ_r spiking **up** in GARD. Given the estimator's ~50× redundancy bias, a plausible mundane | |
| 75 | reading is that a sliding window over rapid coherent compositional drift produces exactly this | |
| 76 | signature — i.e. Φ_r here restates the transition rather than predicting it. Do not over-read the | |
| 77 | sign contrast without an estimator that is symmetric in synergy/redundancy sensitivity. | |
| 78 | ||
| 79 | ## What would actually test the claim | |
| 80 | ||
| 81 | A model with a genuine quiescent phase before a transition. **This was run: see | |
| 82 | [`rung34phi.md`](rung34phi.md).** `rung34.c` does have a lag phase (ignition at gen 13–59 of 500), so | |
| 83 | it hosts the test properly — and Φ_r fails there too, while the raw compositional mean predicts | |
| 84 | ignition timing at p=2.6e-07. The ribozyme null below could be dismissed as a bad test bed; the | |
| 85 | rung34 one cannot. | |
| 86 | ||
| 87 | Still open: `gardenv.c`, the closest structural match to the paper since it *is* GARD, though it is | |
| 88 | one of our nulls (heredity ≈ 0). |
| … | … | @@ -0,0 +1,265 @@ |
| 1 | /* rung34phi.c — Phi_r (causal emergence) probe on rung34.c's compartmentalization emergence. | |
| 2 | * | |
| 3 | * WHY THIS MODEL. The ribozyme probe (ribophi.c) returned a null with a diagnosis: ribozyme.c | |
| 4 | * ignites at generation 6-10, so it has NO quiescent pre-replicative phase and cannot host the | |
| 5 | * paper's "predict the transition from the early trajectory" protocol at all. rung34.c is the | |
| 6 | * better test bed: the compartmentalization trait c starts at 0.10 and climbs to ~0.94 over 500 | |
| 7 | * generations in mutation steps of 0.03, so a genuine lag phase should exist. | |
| 8 | * | |
| 9 | * WHAT IS COPIED VERBATIM from rung34.c: RNG, Dof(), binom(), randn(), and the whole run() loop | |
| 10 | * (within-compartment parasite growth, cooperator->parasite mutation, trait mutation, daughter | |
| 11 | * assortment, survivor selection weighted by cooperator FRACTION). Observation adds no RNG draws. | |
| 12 | * Verified by reproducing rung34.c's single-seed result. | |
| 13 | * | |
| 14 | * SUBSTRATE. Histogram over cooperator FRACTION across the NCOMP compartments (8 equal bins on | |
| 15 | * [0,1]), CLR-transformed with a 0.5 pseudocount, last component dropped -- the analogue of GARD's | |
| 16 | * molecular type counts. Note the trait histogram would be degenerate in the frozen control (c is | |
| 17 | * constant by construction), which is why cooperator fraction is the substrate, not c. | |
| 18 | * | |
| 19 | * PRE-REGISTERED (written before the first run): | |
| 20 | * Q1 WITHIN-ARM (this is the paper's actual design: one parameter set, many seeds, predict | |
| 21 | * timing). In the c-evolvable arm, Phi_r computed on a strictly PRE-IGNITION window correlates | |
| 22 | * across seeds with the generation at which the trait climb ignites. | |
| 23 | * Q2 SPIKES. Igniting runs show Phi_r excursions (>3 sd above the run's own mean) before ignition. | |
| 24 | * Q3 CONTROL -- the one that can kill it. Raw mean cooperator fraction over the SAME window does | |
| 25 | * NOT predict ignition generation. If it does, Phi_r adds nothing. | |
| 26 | * Q4 BETWEEN-ARM. Early Phi_r separates c-evolvable from the frozen control. (Weaker test: the | |
| 27 | * arms differ mechanically from generation 0, so a difference here is not decisive.) | |
| 28 | * | |
| 29 | * Any of these may fail. A null is a result. Do not tune to force a positive. | |
| 30 | * | |
| 31 | * build: cc -O2 -Wall -std=c11 -o rung34phi rung34phi.c -lm | |
| 32 | */ | |
| 33 | #include <math.h> | |
| 34 | #include <stdio.h> | |
| 35 | #include <stdlib.h> | |
| 36 | #include "phir.h" | |
| 37 | ||
| 38 | /* ---------- model constants: verbatim from rung34.c ---------- */ | |
| 39 | #define NCOMP 600 | |
| 40 | #define GEN 500 | |
| 41 | #define SADV 0.6 | |
| 42 | #define DMIN 2 | |
| 43 | #define DMAX 40 | |
| 44 | #define MUTC 0.03 | |
| 45 | #define MUP 0.10 | |
| 46 | ||
| 47 | /* ---------- probe constants ---------- */ | |
| 48 | #define NB 8 /* cooperator-fraction bins */ | |
| 49 | #define P (NB-1) | |
| 50 | #define NSEED 40 | |
| 51 | #define STRIDE 5 | |
| 52 | ||
| 53 | static unsigned long long RNG; | |
| 54 | static unsigned long long xs(void){ RNG^=RNG<<13; RNG^=RNG>>7; RNG^=RNG<<17; return RNG; } | |
| 55 | static double urand(void){ return (double)(xs()>>11)*(1.0/9007199254740992.0); } | |
| 56 | 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); } | |
| 57 | static int Dof(double c){ int D=(int)(DMIN+(1.0-c)*(DMAX-DMIN)+0.5); return D<DMIN?DMIN:(D>DMAX?DMAX:D); } | |
| 58 | static int binom(int n,double p){ int k=0; for(int i=0;i<n;i++) if(urand()<p) k++; return k; } | |
| 59 | ||
| 60 | static double cc[NCOMP]; static int nC[NCOMP]; | |
| 61 | static double dcc[2*NCOMP]; static int dnC[2*NCOMP], dD[2*NCOMP]; static double cum[2*NCOMP]; | |
| 62 | ||
| 63 | /* verbatim rung34 dynamics + observation */ | |
| 64 | static double run(int evolve_c,double *coopf,double cser[GEN],double fser[GEN],double hist[GEN][NB],double histc[GEN][NB]) | |
| 65 | { | |
| 66 | for(int i=0;i<NCOMP;i++){ cc[i]=0.10; nC[i]=Dof(0.10)/2; } | |
| 67 | for(int g=0;g<GEN;g++){ | |
| 68 | for(int i=0;i<NCOMP;i++){ | |
| 69 | int D=Dof(cc[i]); int nc=nC[i]; if(nc>D)nc=D; int np=D-nc; | |
| 70 | double wc=nc, wp=np*(1.0+SADV); double p=(wc+wp>0)?wc/(wc+wp):0; | |
| 71 | int g2=binom(2*D,p); g2=binom(g2,1.0-MUP); double frac=(double)g2/(2*D); | |
| 72 | for(int d=0;d<2;d++){ | |
| 73 | double c2 = evolve_c ? cc[i]+MUTC*randn() : cc[i]; | |
| 74 | if(c2<0)c2=0; | |
| 75 | if(c2>1)c2=1; | |
| 76 | int D2=Dof(c2); | |
| 77 | dcc[2*i+d]=c2; dD[2*i+d]=D2; dnC[2*i+d]=binom(D2,frac); | |
| 78 | } | |
| 79 | } | |
| 80 | double tot=0; for(int i=0;i<2*NCOMP;i++){ tot+=(double)dnC[i]/dD[i]; cum[i]=tot; } | |
| 81 | if(tot<=0){ for(int i=0;i<NCOMP;i++){ nC[i]=0; } } | |
| 82 | else for(int i=0;i<NCOMP;i++){ | |
| 83 | double r=urand()*tot; int lo=0,hi=2*NCOMP-1; | |
| 84 | while(lo<hi){ int mid=(lo+hi)/2; if(cum[mid]<r) lo=mid+1; else hi=mid; } | |
| 85 | cc[i]=dcc[lo]; nC[i]=dnC[lo]; | |
| 86 | } | |
| 87 | /* --- observation only, no RNG consumed --- */ | |
| 88 | { | |
| 89 | double mc=0; long sc=0,sd=0; | |
| 90 | for(int b=0;b<NB;b++) hist[g][b]=0.0; | |
| 91 | for(int i=0;i<NCOMP;i++){ | |
| 92 | int D=Dof(cc[i]); mc+=cc[i]; sc+=nC[i]; sd+=D; | |
| 93 | double fr=(D>0)?(double)nC[i]/D:0.0; | |
| 94 | int b=(int)(fr*NB); if(b>=NB)b=NB-1; if(b<0)b=0; | |
| 95 | hist[g][b]+=1.0; | |
| 96 | } | |
| 97 | if(cser) cser[g]=mc/NCOMP; | |
| 98 | if(fser) fser[g]=(double)sc/sd; | |
| 99 | } | |
| 100 | /* second substrate: distribution of the TRAIT c itself (degenerate in the frozen arm, | |
| 101 | so usable only for the within-arm test, where it is the variable that actually moves) */ | |
| 102 | if(histc){ | |
| 103 | for(int b=0;b<NB;b++) histc[g][b]=0.0; | |
| 104 | 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; } | |
| 105 | } | |
| 106 | } | |
| 107 | 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]); } | |
| 108 | *coopf=(double)sc/sd; return mc/NCOMP; | |
| 109 | } | |
| 110 | ||
| 111 | static double hist_all[2][NSEED][GEN][NB]; | |
| 112 | static double cAll[2][NSEED][GEN], fAll[2][NSEED][GEN]; | |
| 113 | static double clrAll[2][NSEED][GEN][P]; | |
| 114 | static double histc_all[2][NSEED][GEN][NB]; | |
| 115 | static double clrC[2][NSEED][GEN][P]; | |
| 116 | static int ignGen[2][NSEED]; | |
| 117 | static double finalC[2][NSEED], finalF[2][NSEED]; | |
| 118 | ||
| 119 | int main(void) | |
| 120 | { | |
| 121 | phir_seed(0x51ed270b0f8e1c4dULL); | |
| 122 | ||
| 123 | /* ---------- 0. reproduce rung34.c on its own seed ---------- */ | |
| 124 | printf("=== 0. VERIFY: instrumented model reproduces rung34.c (its own seed) ===\n"); | |
| 125 | { | |
| 126 | static double cs[GEN],fs[GEN]; static double hh[GEN][NB],hc[GEN][NB]; double cf; | |
| 127 | RNG=0x9e3779b97f4a7c15ULL; double c1=run(1,&cf,cs,fs,hh,hc); | |
| 128 | printf(" c EVOLVABLE : final mean c = %.3f (D=%d) cooperation frac = %.3f\n",c1,Dof(c1),cf); | |
| 129 | RNG=0x9e3779b97f4a7c15ULL; double c0=run(0,&cf,cs,fs,hh,hc); | |
| 130 | printf(" c FROZEN : final mean c = %.3f (D=%d) cooperation frac = %.3f\n",c0,Dof(c0),cf); | |
| 131 | printf(" (rung34.c reference: evolvable c ~0.94 with cooperation rescued; frozen collapses)\n\n"); | |
| 132 | } | |
| 133 | ||
| 134 | /* ---------- 1. run both arms, NSEED seeds ---------- */ | |
| 135 | const char *armname[2]={"c EVOLVABLE","c FROZEN "}; | |
| 136 | for(int a=0;a<2;a++){ | |
| 137 | for(int s=0;s<NSEED;s++){ | |
| 138 | double cf; | |
| 139 | RNG = 0x9e3779b97f4a7c15ULL ^ (0x100000001b3ULL*(unsigned long long)(s+1)); | |
| 140 | for(int k=0;k<3;k++) xs(); | |
| 141 | finalC[a][s]=run(a==0?1:0,&cf,cAll[a][s],fAll[a][s],hist_all[a][s],histc_all[a][s]); | |
| 142 | finalF[a][s]=cf; | |
| 143 | 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]); } | |
| 144 | int gg=-1; for(int g=0;g<GEN;g++) if(cAll[a][s][g] > 0.10+0.20){ gg=g; break; } | |
| 145 | ignGen[a][s]=gg; | |
| 146 | } | |
| 147 | } | |
| 148 | ||
| 149 | /* ---------- 2. ignition timing: is there a lag phase at all? ---------- */ | |
| 150 | printf("=== 1. Ignition timing (first gen where mean trait c > start+0.20) ===\n"); | |
| 151 | int gmin=GEN; | |
| 152 | for(int a=0;a<2;a++){ | |
| 153 | double gv[NSEED]; int ng=0; int mn=GEN,mx=-1; | |
| 154 | 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]; } | |
| 155 | if(ng){ printf(" %s | ignited %2d/%2d | gen min/median/max = %d / %.1f / %d\n",armname[a],ng,NSEED,mn,median(gv,ng),mx); | |
| 156 | if(mn<gmin) gmin=mn; } | |
| 157 | else printf(" %s | ignited %2d/%2d | (never)\n",armname[a],ng,NSEED); | |
| 158 | printf(" final c = %.3f +- %.3f final coop frac = %.3f +- %.3f\n", | |
| 159 | meanv(finalC[a],NSEED),sdv(finalC[a],NSEED),meanv(finalF[a],NSEED),sdv(finalF[a],NSEED)); | |
| 160 | } | |
| 161 | if(gmin>=GEN){ printf(" no ignition anywhere -- cannot run the test.\n"); return 1; } | |
| 162 | int PRE = gmin; /* gens 0..PRE-1 are pre-ignition in EVERY run */ | |
| 163 | printf(" => earliest ignition anywhere: gen %d. Pre-ignition window = gens 0-%d (%d transitions/run).\n\n", | |
| 164 | gmin,PRE-1,PRE-1); | |
| 165 | ||
| 166 | /* ---------- 3. Q1/Q3: within-arm, does early Phi_r predict ignition timing? ---------- */ | |
| 167 | printf("=== 2. Q1 WITHIN-ARM (c evolvable): does pre-ignition Phi_r predict ignition gen? ===\n"); | |
| 168 | { | |
| 169 | int n=PRE-1; | |
| 170 | if(n<8){ | |
| 171 | printf(" only %d pre-ignition transitions per run -- too few for a per-run estimate.\n",n); | |
| 172 | } else { | |
| 173 | static double pairs[GEN][2*P]; | |
| 174 | double ephi[NSEED], eraw[NSEED], gign[NSEED], fc[NSEED]; int m=0; | |
| 175 | for(int s=0;s<NSEED;s++){ | |
| 176 | if(ignGen[0][s]<0) continue; | |
| 177 | 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]; } | |
| 178 | int ok; double v=phi_r(&pairs[0][0],n,P,&ok); | |
| 179 | if(!ok) continue; | |
| 180 | double q=0; for(int t=0;t<PRE;t++) q+=fAll[0][s][t]; | |
| 181 | ephi[m]=v; eraw[m]=q/PRE; gign[m]=ignGen[0][s]; fc[m]=finalC[0][s]; m++; | |
| 182 | } | |
| 183 | printf(" n=%d runs, window gens 0-%d\n",m,PRE-1); | |
| 184 | printf(" pre-ignition Phi_r : %.5f +- %.5f\n",meanv(ephi,m),sdv(ephi,m)); | |
| 185 | double r1=spearman(ephi,gign,m), r2=spearman(ephi,fc,m); | |
| 186 | double c1=spearman(eraw,gign,m), c2=spearman(eraw,fc,m); | |
| 187 | printf(" Q1 Phi_r vs ignition gen : rho = %+.3f (p=%.3g)\n",r1,spearman_p(r1,m)); | |
| 188 | printf(" Phi_r vs final c : rho = %+.3f (p=%.3g)\n",r2,spearman_p(r2,m)); | |
| 189 | printf(" Q3 raw coopfr vs ignition gen : rho = %+.3f (p=%.3g) <- CONTROL\n",c1,spearman_p(c1,m)); | |
| 190 | printf(" raw coopfr vs final c : rho = %+.3f (p=%.3g) <- CONTROL\n",c2,spearman_p(c2,m)); | |
| 191 | printf(" => Phi_r earns its keep only if |rho| for Phi_r exceeds the control's.\n\n"); | |
| 192 | ||
| 193 | /* Same test on the substrate that actually carries the transition: the TRAIT | |
| 194 | distribution. Cooperator fraction may simply be the wrong variable to watch. */ | |
| 195 | printf(" --- same test, substrate = TRAIT c distribution (the variable that moves) ---\n"); | |
| 196 | double ephiC[NSEED], erawC[NSEED], gignC[NSEED]; int mc2=0; | |
| 197 | for(int s=0;s<NSEED;s++){ | |
| 198 | if(ignGen[0][s]<0) continue; | |
| 199 | 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]; } | |
| 200 | int ok; double v=phi_r(&pairs[0][0],n,P,&ok); | |
| 201 | if(!ok) continue; | |
| 202 | double q=0; for(int t=0;t<PRE;t++) q+=cAll[0][s][t]; | |
| 203 | ephiC[mc2]=v; erawC[mc2]=q/PRE; gignC[mc2]=ignGen[0][s]; mc2++; | |
| 204 | } | |
| 205 | double rC=spearman(ephiC,gignC,mc2), cC=spearman(erawC,gignC,mc2); | |
| 206 | printf(" pre-ignition Phi_r (trait) : %.5f +- %.5f (n=%d)\n",meanv(ephiC,mc2),sdv(ephiC,mc2),mc2); | |
| 207 | printf(" Q1' Phi_r(trait) vs ignition gen : rho = %+.3f (p=%.3g)\n",rC,spearman_p(rC,mc2)); | |
| 208 | printf(" Q3' raw mean c vs ignition gen : rho = %+.3f (p=%.3g) <- CONTROL\n\n",cC,spearman_p(cC,mc2)); | |
| 209 | } | |
| 210 | } | |
| 211 | ||
| 212 | /* ---------- 4. Q4: between-arm on the pre-ignition window ---------- */ | |
| 213 | printf("=== 3. Q4 BETWEEN-ARM on the pre-ignition window (gens 0-%d) ===\n",PRE-1); | |
| 214 | { | |
| 215 | int n=PRE-1; | |
| 216 | if(n>=8){ | |
| 217 | static double pairs[GEN][2*P]; | |
| 218 | double ev[NSEED],fz[NSEED]; int ne=0,nf=0; | |
| 219 | double evR[NSEED],fzR[NSEED]; | |
| 220 | for(int a=0;a<2;a++) for(int s=0;s<NSEED;s++){ | |
| 221 | 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]; } | |
| 222 | int ok; double v=phi_r(&pairs[0][0],n,P,&ok); | |
| 223 | double q=0; for(int t=0;t<PRE;t++) q+=fAll[a][s][t]; q/=PRE; | |
| 224 | if(!ok) continue; | |
| 225 | if(a==0){ ev[ne]=v; evR[ne]=q; ne++; } else { fz[nf]=v; fzR[nf]=q; nf++; } | |
| 226 | } | |
| 227 | printf(" Phi_r evolvable %.5f +- %.5f | frozen %.5f +- %.5f | MW p=%.4g\n", | |
| 228 | meanv(ev,ne),sdv(ev,ne),meanv(fz,nf),sdv(fz,nf),mannwhitney_p(ev,ne,fz,nf)); | |
| 229 | printf(" raw coopfr evolvable %.5f +- %.5f | frozen %.5f +- %.5f | MW p=%.4g <- CONTROL\n", | |
| 230 | meanv(evR,ne),sdv(evR,ne),meanv(fzR,nf),sdv(fzR,nf),mannwhitney_p(evR,ne,fzR,nf)); | |
| 231 | } | |
| 232 | } | |
| 233 | ||
| 234 | /* ---------- 5. Q2: do spikes precede ignition? ---------- */ | |
| 235 | printf("\n=== 4. Q2 SPIKES: Phi_r excursions >3 sd above run mean, before ignition ===\n"); | |
| 236 | { | |
| 237 | int WIN=PRE>25?PRE:25; | |
| 238 | printf(" sliding window = %d gens, stride %d\n",WIN,STRIDE); | |
| 239 | for(int a=0;a<2;a++){ | |
| 240 | int any=0,before=0,nr=0; | |
| 241 | for(int s=0;s<NSEED;s++){ | |
| 242 | static double pairs[GEN][2*P]; static double traj[GEN]; static int tc[GEN]; | |
| 243 | int nw=0; | |
| 244 | for(int st=0; st+WIN<=GEN; st+=STRIDE){ | |
| 245 | int n=WIN-1; | |
| 246 | 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]; } | |
| 247 | int ok; double v=phi_r(&pairs[0][0],n,P,&ok); | |
| 248 | if(ok){ traj[nw]=v; tc[nw]=st+WIN/2; nw++; } | |
| 249 | } | |
| 250 | if(nw<4) continue; | |
| 251 | nr++; | |
| 252 | double mu=meanv(traj,nw), sd=sdv(traj,nw); | |
| 253 | for(int i=0;i<nw;i++) if(traj[i]>mu+3.0*sd){ | |
| 254 | any++; | |
| 255 | if(ignGen[a][s]>=0 && tc[i]<ignGen[a][s]) before++; | |
| 256 | break; | |
| 257 | } | |
| 258 | } | |
| 259 | printf(" %s | runs %2d | with a spike %2d | spike BEFORE ignition %2d\n",armname[a],nr,any,before); | |
| 260 | } | |
| 261 | } | |
| 262 | ||
| 263 | printf("\n=== verdict read off Q1-Q4 above ===\n"); | |
| 264 | return 0; | |
| 265 | } |
| … | … | @@ -0,0 +1,74 @@ |
| 1 | # Φ_r probe against `rung34.c` — the test ribozyme couldn't host, and a replicated null | |
| 2 | ||
| 3 | **Code:** `rung34phi.c` (estimator shared via `phir.h`) · **Companion:** [`ribophi.md`](ribophi.md) | |
| 4 | ||
| 5 | ## Why this model | |
| 6 | ||
| 7 | `ribophi.md` returned a null with a diagnosis: `ribozyme.c` ignites at generation 6–10, so it has no | |
| 8 | quiescent pre-transition phase and cannot host the paper's "predict the transition from the early | |
| 9 | trajectory" protocol at all. `rung34.c` should be different — the compartmentalization trait c climbs | |
| 10 | from 0.10 to ~0.94 in mutation steps of 0.03 over 500 generations. | |
| 11 | ||
| 12 | **It is different.** Measured ignition (first gen where mean c exceeds start+0.20): **gen 13–59, | |
| 13 | median 28**, in 40/40 runs, out of 500 generations. A genuine lag phase exists. The test can run. | |
| 14 | ||
| 15 | Model reproduced verbatim: evolvable c → 0.942 (D=4), cooperation 0.873; frozen → 0.100, cooperation | |
| 16 | 0.244. Matches `rung34.c`. | |
| 17 | ||
| 18 | ## Design — this time matching the paper's actual shape | |
| 19 | ||
| 20 | The paper's design is *one parameter set, many seeds, predict the timing of the transition*. So the | |
| 21 | primary test here is **within-arm across seeds**: in the c-evolvable arm, does Φ_r computed on a | |
| 22 | strictly pre-ignition window (gens 0–12, before *any* run has ignited) predict which seeds ignite | |
| 23 | sooner? Two substrates, because the choice matters: | |
| 24 | ||
| 25 | - cooperator-fraction histogram (usable in both arms; the frozen arm's trait histogram is degenerate) | |
| 26 | - **trait-c histogram** — the variable that actually moves, usable within-arm | |
| 27 | ||
| 28 | ## Results | |
| 29 | ||
| 30 | | test | Φ_r | control (raw mean) | | |
| 31 | |---|---|---| | |
| 32 | | Q1 coop-fraction substrate vs ignition gen | ρ = −0.102 (p=0.53) | ρ = −0.134 (p=0.41) | | |
| 33 | | Q1′ **trait substrate vs ignition gen** | ρ = **+0.161 (p=0.32)** | ρ = **−0.641 (p=2.6e-07)** | | |
| 34 | | Q4 between-arm, evolvable vs frozen | p=0.46 | p=0.14 | | |
| 35 | | Q2 spikes before ignition | **0 of 80 runs** | — | | |
| 36 | ||
| 37 | **Q1 fails. Q2 fails. Q4 fails. Q3 (the control) is decisive.** On the substrate that actually carries | |
| 38 | the transition, the raw compositional mean predicts ignition timing at **p=2.6e-07** while Φ_r on the | |
| 39 | same data predicts nothing (p=0.32). Runs whose trait mean drifts up faster in the first 12 | |
| 40 | generations ignite sooner — trivially, mechanistically predictive. Φ_r adds nothing on top. | |
| 41 | ||
| 42 | Note the coop-fraction substrate gives a *double* null: neither Φ_r nor the control predicts. That is | |
| 43 | the right sanity check — it shows the trait-substrate result is not an artifact of the method but a | |
| 44 | genuine signal that Φ_r fails to capture and a simple mean captures easily. | |
| 45 | ||
| 46 | ## Combined verdict across both probes | |
| 47 | ||
| 48 | | | `ribozyme.c` | `rung34.c` | | |
| 49 | |---|---|---| | |
| 50 | | lag phase before transition | **none** (ignites gen 6–10) | **yes** (gen 13–59 of 500) | | |
| 51 | | Φ_r predicts the transition | no | no | | |
| 52 | | raw composition predicts it | yes, p=3.1e-17 | yes, p=2.6e-07 | | |
| 53 | | Φ_r spikes before ignition | 0 / 100 runs | 0 / 80 runs | | |
| 54 | ||
| 55 | Two models, one of which genuinely hosts the test, and in both the simple compositional baseline beats | |
| 56 | Φ_r decisively while Φ_r contributes nothing. The ribozyme null could be dismissed as a bad test bed. | |
| 57 | This one cannot. | |
| 58 | ||
| 59 | ## What this does and does not license | |
| 60 | ||
| 61 | **Does not** show the paper is wrong. Three honest differences remain, any of which could account for it: | |
| 62 | ||
| 63 | 1. **Different models.** Neither of ours is GARD. GARD is accretion-and-fission of a single assembly | |
| 64 | with no selection; ours are population models under selection from generation 1. | |
| 65 | 2. **Estimator.** Ours is Gaussian log-determinant ΦID with a measured ~50× bias toward redundancy | |
| 66 | over synergy (`ribophi.c` section 0a). The paper's preprocessing implies Gaussian too, but their | |
| 67 | ΦID/PID implementation is not published yet ("code available upon publication") and may differ. | |
| 68 | 3. **Protocol.** They fit an MLP on the Φ_r trajectory and beat composition/flux baselines; we use a | |
| 69 | rank correlation against a single scalar mean. Ours is the harder baseline to beat — if a single | |
| 70 | mean beats Φ_r, that is informative — but it is not the same test. | |
| 71 | ||
| 72 | **Does** establish that in two independent non-GARD models with known ground-truth transitions, Φ_r | |
| 73 | as we can best reconstruct it is not a leading indicator, and a trivial baseline is. Worth putting to | |
| 74 | the authors alongside the question about boundaries. |