/* gardphi.c — Phi_r probe on OUR GARD (gardenv.c). The only model of ours that is actually GARD.
*
* WHY THIS ONE MATTERS. ribophi.c and rung34phi.c tested Pigozzi & Levin's causal-emergence claim
* (arXiv 2607.28250) by ANALOGY, in models that are not GARD. Both returned nulls, and the honest
* caveat in both was the same: different model, so the disagreement could sit anywhere. gardenv.c
* IS the GARD family -- token types, a catalysis matrix beta, accretion, split-and-regrow,
* compositional (not sequence) inheritance. So this probe can speak to their result directly.
*
* MATCHING THEIR SETUP, not ours:
* - SINGLE LINEAGE. The paper traces one assembly, fissioning and continuing from one daughter.
* gardenv.c normally runs 200 pockets with between-pocket selection; here selection is OFF and
* one lineage is followed, because the paper is explicit that GARD "is essentially a stochastic
* model of events before natural selection takes place".
* - MOLECULAR TIME, not generations. They compute Phi_r over molecular steps. ribophi.c sampled
* once per generation and that coarseness was part of why it could not host the test. Here the
* composition is sampled every SAMPLE join/leave events, giving ~9000 points per run.
* - THEIR SELF-REPLICATOR CRITERION. "Recurring compositions inherited across generations": find
* the most recurring composition, then a timestep is in-self-replicator if its cosine to that
* composition exceeds a threshold. A run enters and leaves self-replication over its lifetime.
*
* DIMENSION REDUCTION (declared, because it is a real deviation). They use Ng=100 parts and an
* unspecified minimum-information-bipartition search; exhaustive search over 100 parts is
* impossible, so they must use a heuristic they do not describe. We reduce N=40 token types to 8
* coordinates -- the 7 most abundant types plus an "other" bin -- chosen from the FIRST 25% of each
* run only, so no information from the predicted portion leaks into the coordinate choice. Then
* CLR, drop one, 7 parts, exhaustive bipartition (63) as in the other probes.
*
* PRE-REGISTERED:
* R1 (their Fig 3/4, the headline): Phi_r is HIGHER during timesteps in self-replication than
* outside it. Per run by Mann-Whitney; report how many runs are significantly positive.
* R2 (their Fig 5): mean Phi_r over the FIRST 25% of a run predicts the fraction of the REMAINING
* 75% spent in self-replication, across runs (Spearman).
* R3 CONTROL for R2, the one that can kill it: the same prediction from a raw compositional
* statistic over the same first 25% (mean drift = 1 - cosine between consecutive samples).
* Phi_r earns its keep only if it beats this. NOTE R1 is not circular -- Phi_r is not derived
* from the label -- but R2's control must be, since the label itself is compositional.
* R4 TEST OF THE TEST: does self-replication actually occur and fluctuate? If a run is in-replicator
* ~always or ~never, there is nothing to predict and that run is uninformative. Report the
* distribution of in-replicator fraction before reading anything else.
*
* Prior from our own arc: gardenv.c measured heredity ~0 in this chemistry -- composomes do not
* breed true. So compositions may satisfy the paper's similarity criterion while transmitting
* nothing. If R1 holds while heredity is ~0, that is not a contradiction of the paper; it locates
* precisely what their "self-replication" is tracking.
*
* build: cc -O2 -Wall -std=c11 -o gardphi gardphi.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "phir.h"
/* ---------- verbatim from gardenv.c ---------- */
#define N 40
#define S 150
#define GEN 60
#define MAT 1500
#define PCAT 0.15
#define RHO0 0.4
/* ---------- probe ---------- */
#define SAMPLE 10 /* sample the composition every SAMPLE molecular events */
#define MAXT 12000
#define NB 8
#define P (NB-1)
#define WINL 200 /* Phi_r window, in samples */
#define STRIDE 50
#define NRUN 30
#define SIMTHR 0.90 /* cosine to the most recurring composition = "in self-replicator" */
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 randnv(void){ double u=urand(),v=urand(); if(u<1e-300)u=1e-300; return sqrt(-2*log(u))*cos(6.283185307*v); }
static double bta[N][N], rate[N];
static void build_chem(void){ for(int i=0;i<N;i++)for(int j=0;j<N;j++) bta[i][j]=(i!=j&&urand()<PCAT)? exp(1.8*randnv())*0.30 : 0.0; }
static void recompute(const int*n){ for(int i=0;i<N;i++){ double r=RHO0; for(int j=0;j<N;j++) r+=bta[i][j]*n[j]; rate[i]=r; } }
static int pick_add(void){ double t=0; for(int i=0;i<N;i++)t+=rate[i]; double u=urand()*t,a=0; for(int i=0;i<N;i++){a+=rate[i];if(u<=a)return i;} return N-1; }
static int pick_rem(const int*n,int sz){ double u=urand()*sz,a=0; for(int j=0;j<N;j++){a+=n[j];if(u<=a)return j;} return N-1; }
static void addm(int*n,int i){ n[i]++; for(int k=0;k<N;k++) rate[k]+=bta[k][i]; }
static void remm(int*n,int j){ n[j]--; for(int k=0;k<N;k++) rate[k]-=bta[k][j]; }
static double cosv(const double*a,const double*b){ double d=0,na=0,nb=0; for(int i=0;i<N;i++){ d+=a[i]*b[i]; na+=a[i]*a[i]; nb+=b[i]*b[i]; } if(na<=0||nb<=0)return 0; return d/(sqrt(na)*sqrt(nb)); }
static double traj[MAXT][N];
static int ntraj;
static void snap(const int*n){ if(ntraj>=MAXT) return; for(int i=0;i<N;i++) traj[ntraj][i]=n[i]; ntraj++; }
/* mature + grow, sampling composition in MOLECULAR time (the paper's substrate) */
static void mature_s(int*n,int sz){ recompute(n); for(int t=0;t<MAT;t++){ int j=pick_rem(n,sz); remm(n,j); addm(n,pick_add()); if(t%SAMPLE==0) snap(n); } }
static void grow_to_s(int*n,int*sz,int tgt){ recompute(n); int c=0; while(*sz<tgt){ addm(n,pick_add()); (*sz)++; if((c++)%SAMPLE==0) snap(n); } }
/* one lineage: accrete, mature, fission, continue from a daughter. NO selection. */
static void run_lineage(void)
{
int n[N]; memset(n,0,sizeof n);
for(int s=0;s<S;s++) n[(int)(xs()%N)]++;
ntraj=0;
mature_s(n,S);
for(int g=0;g<GEN;g++){
int D[N]; memset(D,0,sizeof D); int ds=0;
for(int i=0;i<N;i++) for(int c=0;c<n[i];c++) if(urand()<0.5){ D[i]++; ds++; }
if(ds<1){ D[0]=1; ds=1; }
grow_to_s(D,&ds,S);
mature_s(D,S);
memcpy(n,D,sizeof D);
}
}
int main(void)
{
phir_seed(0x9e3779b97f4a7c15ULL);
printf("=== Phi_r on OUR GARD (gardenv.c chemistry): single lineage, no selection ===\n");
printf(" %d runs, %d tokens, %d generations, composition sampled every %d molecular events.\n",NRUN,N,GEN,SAMPLE);
printf(" Self-replicator = cosine to the most recurring composition > %.2f (their criterion).\n",SIMTHR);
printf(" Coordinates: 7 most abundant token types + other, chosen from the FIRST 25%% only.\n\n");
double earlyPhi[NRUN], lateFrac[NRUN], earlyDrift[NRUN], inFrac[NRUN];
int nrun=0, nsig_pos=0, nsig_neg=0, ntested=0;
for(int r=0;r<NRUN;r++){
RNG = 0x9e3779b97f4a7c15ULL ^ (0x100000001b3ULL*(unsigned long long)(r+1));
for(int k=0;k<5;k++) xs();
build_chem();
run_lineage();
if(ntraj < 4*WINL) continue;
int q1 = ntraj/4; /* first 25% */
/* --- coordinates from the first 25% only: 7 most abundant types + other --- */
double ab[N]; for(int i=0;i<N;i++) ab[i]=0;
for(int t=0;t<q1;t++) for(int i=0;i<N;i++) ab[i]+=traj[t][i];
int top[NB-1];
{ int used[N]; for(int i=0;i<N;i++) used[i]=0;
for(int s=0;s<NB-1;s++){ int bi=-1; double bv=-1;
for(int i=0;i<N;i++) if(!used[i] && ab[i]>bv){ bv=ab[i]; bi=i; }
top[s]=bi; used[bi]=1; } }
static double clr[MAXT][P];
for(int t=0;t<ntraj;t++){
double bins[NB]; for(int b=0;b<NB;b++) bins[b]=0;
double tot=0; for(int i=0;i<N;i++) tot+=traj[t][i];
for(int s=0;s<NB-1;s++) bins[s]=traj[t][top[s]];
double used=0; for(int s=0;s<NB-1;s++) used+=bins[s];
bins[NB-1]=tot-used;
phir_clr(bins,NB,clr[t]);
}
/* --- most recurring composition: biggest cosine>0.85 cluster over the whole run --- */
double cent[N]; int bestcount=-1;
{
int NPROBE=60;
for(int a=0;a<NPROBE;a++){
int idx=(int)((double)a/NPROBE*ntraj);
int cnt=0; for(int t=0;t<ntraj;t+=5) if(cosv(traj[idx],traj[t])>0.85) cnt++;
if(cnt>bestcount){ bestcount=cnt; memcpy(cent,traj[idx],sizeof cent); }
}
}
static int lab[MAXT];
int nin=0;
for(int t=0;t<ntraj;t++){ lab[t] = (cosv(traj[t],cent)>SIMTHR); nin+=lab[t]; }
inFrac[nrun]=(double)nin/ntraj;
/* --- Phi_r trajectory in molecular time --- */
static double pairs[MAXT][2*P];
static double pv[MAXT]; static int pc[MAXT]; int nw=0;
for(int st=0; st+WINL<=ntraj; st+=STRIDE){
int n2=WINL-1;
for(int t=0;t<n2;t++) for(int i=0;i<P;i++){ pairs[t][i]=clr[st+t][i]; pairs[t][P+i]=clr[st+t+1][i]; }
int ok; double v=phi_r(&pairs[0][0],n2,P,&ok);
if(ok){ pv[nw]=v; pc[nw]=st+WINL/2; nw++; }
}
if(nw<8) continue;
/* --- R1: Phi_r inside vs outside self-replication (per run) --- */
double insv[MAXT], outv[MAXT]; int ni=0,no=0;
for(int i=0;i<nw;i++){ if(lab[pc[i]]) insv[ni++]=pv[i]; else outv[no++]=pv[i]; }
if(ni>=4 && no>=4){
ntested++;
double p=mannwhitney_p(insv,ni,outv,no);
double mi=meanv(insv,ni), mo=meanv(outv,no);
if(p<0.05){ if(mi>mo) nsig_pos++; else nsig_neg++; }
}
/* --- R2/R3: early Phi_r and early drift vs later self-replication --- */
{
double sPhi=0; int cPhi=0;
for(int i=0;i<nw;i++) if(pc[i]<q1){ sPhi+=pv[i]; cPhi++; }
earlyPhi[nrun] = cPhi? sPhi/cPhi : NAN;
double d=0; int cd=0;
for(int t=1;t<q1;t++){ d += 1.0-cosv(traj[t-1],traj[t]); cd++; }
earlyDrift[nrun] = cd? d/cd : NAN;
int late=0,lc=0;
for(int t=q1;t<ntraj;t++){ late+=lab[t]; lc++; }
lateFrac[nrun] = lc? (double)late/lc : NAN;
}
nrun++;
}
/* ---------- R4 first: did self-replication happen and fluctuate? ---------- */
printf(" R4 TEST OF THE TEST — in-self-replicator fraction across %d runs:\n",nrun);
{
double f[NRUN]; int nf=0, degenerate=0;
for(int i=0;i<nrun;i++){ f[nf++]=inFrac[i]; if(inFrac[i]<0.02||inFrac[i]>0.98) degenerate++; }
printf(" mean %.3f, median %.3f, min %.3f, max %.3f | degenerate runs (always/never): %d/%d\n",
meanv(f,nf), median(f,nf), f[0], f[nf-1], degenerate, nrun);
double mn=1,mx=0; for(int i=0;i<nf;i++){ if(f[i]<mn)mn=f[i]; if(f[i]>mx)mx=f[i]; }
printf(" range %.3f .. %.3f\n",mn,mx);
if(degenerate>nrun/2)
printf(" => MOSTLY DEGENERATE: little to predict. Treat everything below as uninformative.\n");
else
printf(" => self-replication occurs and fluctuates; the test is live.\n");
}
/* ---------- R1 ---------- */
printf("\n R1 (their headline): Phi_r higher DURING self-replication than outside?\n");
printf(" runs tested %d | significantly HIGHER %d | significantly LOWER %d | n.s. %d\n",
ntested,nsig_pos,nsig_neg,ntested-nsig_pos-nsig_neg);
/* ---------- R2 / R3 ---------- */
printf("\n R2/R3: early (first 25%%) predictor vs later self-replication fraction, across runs\n");
{
double a[NRUN],b[NRUN],c[NRUN]; int m=0;
for(int i=0;i<nrun;i++){
if(isnan(earlyPhi[i])||isnan(lateFrac[i])||isnan(earlyDrift[i])) continue;
a[m]=earlyPhi[i]; b[m]=lateFrac[i]; c[m]=earlyDrift[i]; m++;
}
double r1=spearman(a,b,m), r2=spearman(c,b,m);
printf(" n=%d runs\n",m);
printf(" R2 early Phi_r vs later self-replication : rho = %+.3f (p=%.3g)\n",r1,spearman_p(r1,m));
printf(" R3 early drift vs later self-replication : rho = %+.3f (p=%.3g) <- CONTROL\n",r2,spearman_p(r2,m));
printf(" => Phi_r earns its keep only if |rho| exceeds the control's.\n");
}
printf("\n=== verdict read off R1-R4 ===\n");
return 0;
}