/* rafcomp.c — TEST B: chain rung 1 into rung 4. Does compartment-level selection still work when
* catalytic quality is a COLLECTIVE closure property instead of a per-unit one?
*
* WHY. SYNTHESIS's standing admission is that no rung ever produced the next — every rung is an
* isolated island. chemoton.c derives compartment-level selection from growth coupling, but it still
* SUPPLIES phi as a per-unit fold->catalysis map. rung1.c separately shows autocatalytic sets (RAFs)
* self-organising from random chemistry. This chains them: catalytic output is read off an EMERGENT
* RAF rather than a supplied fold rule, putting rung 1 and rung 4 in ONE model.
*
* THE REASON IT MIGHT FAIL, AND WHY THAT IS THE POINT. phi is intrinsic to a unit: a compartment's
* quality is the average of its units' qualities, so assortment shuffles good units around without
* destroying what makes them good. RAF membership is NOT intrinsic — it is a property of the SET.
* A reaction is in the RAF only if its reactants are producible and it is catalyzed by something
* producible, all within the set present. So splitting a compartment can DESTROY the closure that
* made the parent productive, even though every individual reaction survives.
*
* That collides head-on with Test A, which measured that ASSORTMENT VARIANCE is what purges
* parasites (propagule beats fission by +1.14 phi at matched realized size). Variance is needed to
* purge parasites; variance destroys collective closure. Test B asks which wins.
*
* WHAT IS SUPPLIED vs WHAT MUST EMERGE:
* SUPPLIED: a random catalytic reaction network (rung1.c's given chemistry), a food set, the
* growth/division/death machinery from chemoton.c, and a conserved material pool. Catalysis
* density is SWEPT.
* MUST EMERGE: the RAF itself (nothing designates which reactions are "good" — closure is computed,
* not assigned), and compartment-level differential reproduction (nothing samples parents by
* fitness; compartments divide because they grew).
* UNTOUCHED: the tragedy. Reactions OUTSIDE the RAF are cheap — they carry no catalytic obligation
* — so they replicate faster, exactly as low-phi parasites did.
*
* PRE-REGISTERED:
* T1 TEST OF THE TEST: do RAFs form inside compartments at all? Report mean maxRAF size against the
* random-subset baseline. If no compartment ever closes a RAF, nothing was selected on and a
* null says nothing.
* T2 coupled uptake raises mean RAF size above the random baseline. Compartment selection survives
* the move from per-unit to collective quality.
* T3 CONTROL: uptake blind to RAF size, everything else identical. Must stay at baseline.
* T4 THE REVERSAL: propagule (high variance) should be WORSE than fission here, opposite to Test A,
* because a single-unit bottleneck cannot carry a closure. If propagule still wins, collective
* function is more robust to assortment than expected and Test A's finding generalises.
*
* Do not tune to force a positive. A null with T1 passing is a real result: it would say the ladder
* does not chain, and that chemoton's success depended on quality being per-unit.
*
* build: cc -O2 -Wall -std=c11 -o rafcomp rafcomp.c -lm
*/
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
/* ---- reaction network (rung1.c shape, scaled so per-compartment RAF is cheap to compute) ---- */
#define NMOL 30
#define NFOOD 6
#define MRXN 60
/* ---- compartment machinery (chemoton.c shape) ---- */
#define MAXCOMP 2400
#define MAXSZ 24
#define GEN 200
#define BETA 2.0
#define COSTUNIT 1.0
#define TOTALMAT 4000.0
#define UPTAKE 0.55
#define DEATH 0.10
#define MUTR 0.02 /* chance a copied reaction is replaced by a random one */
static unsigned long long RNG;
static unsigned long long xs(void){ RNG^=RNG<<13; RNG^=RNG>>7; RNG^=RNG<<17; return RNG; }
static double urand(void){ return (double)(xs()>>11)*(1.0/9007199254740992.0); }
static int rint_(int n){ return (int)(xs()%(unsigned)n); }
static int r1[MRXN],r2[MRXN],prod[MRXN];
static unsigned char cat[MRXN][NMOL];
static void build_net(double p){
for(int j=0;j<MRXN;j++){
r1[j]=rint_(NMOL); r2[j]=rint_(NMOL);
prod[j]=NFOOD+rint_(NMOL-NFOOD);
for(int m=0;m<NMOL;m++) cat[j][m]=(urand()<p)?1:0;
}
}
/* maxRAF over a SUBSET of reactions (prune-to-fixpoint, rung1.c's algorithm).
inraf[] is filled with 1 for reactions surviving in the closure. Returns its size. */
static int max_raf_subset(const int *rx,int n,unsigned char *inraf)
{
unsigned char act[MAXSZ], W[NMOL];
for(int i=0;i<n;i++) act[i]=1;
for(;;){
for(int m=0;m<NMOL;m++) W[m]=(m<NFOOD)?1:0;
int grew=1;
while(grew){ grew=0;
for(int i=0;i<n;i++){ if(!act[i]) continue; int j=rx[i];
if(!W[r1[j]]||!W[r2[j]]||W[prod[j]]) continue;
int cz=0; for(int m=0;m<NMOL;m++) if(cat[j][m]&&W[m]){cz=1;break;}
if(cz){ W[prod[j]]=1; grew=1; }
}
}
int changed=0;
for(int i=0;i<n;i++){ if(!act[i]) continue; int j=rx[i];
int sup=W[r1[j]]&&W[r2[j]];
int cz=0; for(int m=0;m<NMOL;m++) if(cat[j][m]&&W[m]){cz=1;break;}
if(!(sup&&cz)){ act[i]=0; changed=1; }
}
if(!changed) break;
}
int c=0; for(int i=0;i<n;i++){ inraf[i]=act[i]; c+=act[i]; }
return c;
}
typedef struct { int rx[MAXSZ]; unsigned char inraf[MAXSZ]; int n; int raf; double mat; } Comp;
static Comp C[MAXCOMP];
static int ncomp;
static double pool;
static double demand[MAXCOMP];
static double run(int coupled,int divsz,int propagule,double *raf_end,double *ncomp_end,
double *msize,double *ndiv,double *matcheck)
{
int seed_n = divsz/2>2?divsz/2:2;
ncomp=600;
for(int c=0;c<ncomp;c++){
C[c].n=seed_n; C[c].mat=0;
for(int k=0;k<seed_n;k++) C[c].rx[k]=rint_(MRXN);
}
{ long u=0; for(int c=0;c<ncomp;c++) u+=C[c].n; pool=TOTALMAT-u*COSTUNIT; }
double totdiv=0;
for(int g=0;g<GEN;g++){
for(int c=0;c<ncomp;c++) C[c].raf = max_raf_subset(C[c].rx,C[c].n,C[c].inraf);
double D=0;
for(int c=0;c<ncomp;c++){
double out = coupled ? (double)C[c].raf : 0.5*C[c].n; /* control: blind to closure */
demand[c]=UPTAKE*out; D+=demand[c];
}
double frac=(D>0&&D>pool)?pool/D:1.0;
for(int c=0;c<ncomp;c++){ double t=demand[c]*frac; C[c].mat+=t; pool-=t; }
/* growth: reactions OUTSIDE the RAF are cheap and replicate faster (the tragedy) */
for(int c=0;c<ncomp;c++){
while(C[c].mat>=COSTUNIT && C[c].n<MAXSZ && C[c].n>0){
double tot=0,cum[MAXSZ];
for(int k=0;k<C[c].n;k++){ tot += 1.0 + BETA*(C[c].inraf[k]?0.0:1.0); cum[k]=tot; }
double r=urand()*tot; int lo=0,hi=C[c].n-1;
while(lo<hi){ int mid=(lo+hi)/2; if(cum[mid]<r) lo=mid+1; else hi=mid; }
int nj = (urand()<MUTR)? rint_(MRXN) : C[c].rx[lo];
C[c].rx[C[c].n]=nj; C[c].inraf[C[c].n]=0; C[c].n++; C[c].mat-=COSTUNIT;
}
}
for(int c=0;c<ncomp;c++)
for(int k=0;k<C[c].n;){
if(urand()<DEATH){ C[c].rx[k]=C[c].rx[C[c].n-1]; C[c].inraf[k]=C[c].inraf[C[c].n-1]; C[c].n--; pool+=COSTUNIT; }
else k++;
}
int base=ncomp;
for(int c=0;c<base;c++){
if(C[c].n>=divsz && ncomp<MAXCOMP){
Comp *p=&C[c],*q=&C[ncomp];
if(propagule){
int pick=rint_(p->n); double sh=p->mat/(double)p->n;
q->n=0; q->rx[q->n]=p->rx[pick]; q->inraf[q->n]=0; q->n++; q->mat=sh; p->mat-=sh;
p->rx[pick]=p->rx[p->n-1]; p->inraf[pick]=p->inraf[p->n-1]; p->n--;
} else {
q->n=0; q->mat=p->mat*0.5; p->mat*=0.5; int w=0;
for(int k=0;k<p->n;k++){
if(urand()<0.5){ q->rx[q->n]=p->rx[k]; q->inraf[q->n]=p->inraf[k]; q->n++; }
else { p->rx[w]=p->rx[k]; p->inraf[w]=p->inraf[k]; w++; }
}
p->n=w;
}
ncomp++; totdiv++;
}
}
for(int c=0;c<ncomp;){ if(C[c].n<=0){ pool+=C[c].mat; C[c]=C[ncomp-1]; ncomp--; } else c++; }
if(ncomp<=0) break;
}
double sr=0; long u=0; double stored=0;
for(int c=0;c<ncomp;c++){
C[c].raf=max_raf_subset(C[c].rx,C[c].n,C[c].inraf);
sr+=C[c].raf; u+=C[c].n; stored+=C[c].mat;
}
*raf_end = ncomp? sr/ncomp : 0;
*ncomp_end=ncomp; *msize= ncomp? (double)u/ncomp : 0;
*ndiv=totdiv/GEN; *matcheck=pool+stored+u*COSTUNIT;
return *raf_end;
}
/* baseline: expected maxRAF of a RANDOM subset of the same size, same network */
static double random_baseline(int n,int reps)
{
double s=0; int rx[MAXSZ]; unsigned char ir[MAXSZ];
for(int t=0;t<reps;t++){ for(int k=0;k<n;k++) rx[k]=rint_(MRXN); s+=max_raf_subset(rx,n,ir); }
return s/reps;
}
int main(void)
{
unsigned long long seeds[5];
seeds[0]=0x9e3779b97f4a7c15ULL; seeds[1]=0xd1b54a32d192ed03ULL; seeds[2]=0x2545f4914f6cdd1dULL;
for(int i=3;i<5;i++) seeds[i]=0x9e3779b97f4a7c15ULL ^ (0x100000001b3ULL*(unsigned long long)(i+1));
int NS=5;
printf("=== TEST B: rung 1 chained into rung 4 — collective quality instead of per-unit ===\n");
printf(" Catalytic output = size of the maxRAF of the reactions a compartment actually holds.\n");
printf(" Nothing designates which reactions are good: closure is COMPUTED, not assigned.\n");
printf(" Reactions outside the RAF carry no catalytic obligation, so they replicate %.0fx faster.\n",1+BETA);
printf(" THE TENSION: Test A found assortment VARIANCE purges parasites (+1.14 at matched size).\n");
printf(" But a RAF is a property of the SET, so variance should DESTROY closure. Which wins?\n\n");
/* ---- MATCHED REALIZED SIZE. Test B's reversal was established at matched division
THRESHOLD, but propagule mode also shrinks realized compartment size, and in this model
bigger compartments hold bigger RAFs. So fission > propagule may be a pure size artifact
with no closure destruction at all. Sweep division size finely at one catalysis density,
record each arm's OWN realized size, and interpolate fission onto propagule's sizes. ---- */
double dens = 0.10;
int sizes[]={4,6,8,10,12,16,20,24};
int NZ=(int)(sizeof(sizes)/sizeof(sizes[0]));
double fr[16], fsz[16], pr[16], psz[16], cf[16], cp[16], base[16];
int any_raf=0, ctrl_flat=1;
printf(" catalysis density %.2f, %d seeds\n\n", dens, NS);
printf(" divsz | rand base | FISSION raf / size | PROPAGULE raf / size | fis ctrl | prop ctrl\n");
printf(" ------+-----------+--------------------+----------------------+----------+----------\n");
for(int z=0;z<NZ;z++){
double bf=0,bs=0,pf=0,ps=0,c1=0,c2=0,bb=0;
for(int s=0;s<NS;s++){
double re,n,m,dv,mm;
RNG=seeds[s]; build_net(dens); bb+=random_baseline(sizes[z]/2>2?sizes[z]/2:2,200);
RNG=seeds[s]; build_net(dens); bf+=run(1,sizes[z],0,&re,&n,&m,&dv,&mm); bs+=m;
RNG=seeds[s]; build_net(dens); pf+=run(1,sizes[z],1,&re,&n,&m,&dv,&mm); ps+=m;
RNG=seeds[s]; build_net(dens); c1+=run(0,sizes[z],0,&re,&n,&m,&dv,&mm);
RNG=seeds[s]; build_net(dens); c2+=run(0,sizes[z],1,&re,&n,&m,&dv,&mm);
}
fr[z]=bf/NS; fsz[z]=bs/NS; pr[z]=pf/NS; psz[z]=ps/NS; cf[z]=c1/NS; cp[z]=c2/NS; base[z]=bb/NS;
printf(" %2d | %5.2f | %5.2f / %5.2f | %5.2f / %5.2f | %5.2f | %5.2f\n",
sizes[z],base[z],fr[z],fsz[z],pr[z],psz[z],cf[z],cp[z]);
if(fr[z]>0.2||pr[z]>0.2) any_raf=1;
if(cf[z]>base[z]+0.3||cp[z]>base[z]+0.3) ctrl_flat=0;
}
/* fission's realized size is NON-MONOTONIC in division size (it rises then falls as the arm
collapses at large divsz), so the fission points must be SORTED BY SIZE before interpolating.
The first version bounds-checked against fsz[last] and silently dropped almost every point. */
printf("\n Does RAF track SIZE regardless of reproduction mode? (per-cell, same divsz)\n");
int tracks=0;
for(int z=0;z<NZ;z++){
int bigger_f = fsz[z]>psz[z];
int rafbig_f = fr[z]>pr[z];
if(bigger_f==rafbig_f) tracks++;
printf(" divsz %2d : larger arm = %s, larger RAF = %s %s\n", sizes[z],
bigger_f?"fission ":"propagule", rafbig_f?"fission ":"propagule",
(bigger_f==rafbig_f)?"tracks":"DIVERGES");
}
printf(" RAF tracks realized size in %d/%d cells\n",tracks,NZ);
printf("\n MATCHED REALIZED SIZE — propagule vs fission interpolated (fission sorted by size):\n");
int ord[16]; for(int z=0;z<NZ;z++) ord[z]=z;
for(int a=1;a<NZ;a++){ int k=ord[a],b=a-1; while(b>=0&&fsz[ord[b]]>fsz[k]){ ord[b+1]=ord[b]; b--; } ord[b+1]=k; }
double sumd=0; int nd=0;
for(int z=0;z<NZ;z++){
double t=psz[z];
if(t<=fsz[ord[0]]||t>=fsz[ord[NZ-1]]) continue;
int lo=-1; for(int y=0;y+1<NZ;y++) if(fsz[ord[y]]<=t&&fsz[ord[y+1]]>=t){ lo=y; break; }
if(lo<0) continue;
double d0=fsz[ord[lo]], d1=fsz[ord[lo+1]];
double w=(d1>d0)?(t-d0)/(d1-d0):0.0;
double fi=fr[ord[lo]]+w*(fr[ord[lo+1]]-fr[ord[lo]]);
double d=pr[z]-fi;
printf(" size %5.2f : propagule %5.2f vs fission-interp %5.2f diff %+.2f\n", t, pr[z], fi, d);
sumd+=d; nd++;
}
double mean_adv = nd? sumd/nd : 0.0/0.0;
printf(" mean propagule advantage at MATCHED size: %+.2f over %d points\n\n", mean_adv, nd);
printf(" T1 RAFs formed : %s\n", any_raf?"YES":"NO — test did not run");
printf(" T3 CONTROL blind to closure stayed flat : %s\n", ctrl_flat?"YES":"NO — CONFOUNDED");
printf(" T4 at matched DIVISION THRESHOLD, fission won in the earlier run.\n");
printf(" T4' at matched REALIZED SIZE, propagule advantage = %+.2f\n\n", mean_adv);
if(!any_raf) printf(" => TEST DID NOT RUN.\n");
else if(!ctrl_flat) printf(" => CONFOUNDED by the control.\n");
else if(mean_adv < -0.2)
printf(" => THE REVERSAL SURVIVES. Even at matched realized size the propagule bottleneck is\n"
" WORSE for collective closure (%+.2f). So it is not a size artifact: a single-unit\n"
" bottleneck genuinely cannot carry a RAF, and Test A's finding does NOT generalise.\n"
" Optimal assortment is a trade-off set by whether quality is per-unit or collective.\n", mean_adv);
else if(mean_adv > 0.2)
printf(" => THE REVERSAL WAS AN ARTIFACT, AND IT INVERTS. At matched realized size the propagule\n"
" bottleneck is BETTER (%+.2f), same direction as Test A. The earlier fission win was\n"
" entirely propagule mode shrinking compartments, and in this model smaller compartments\n"
" hold smaller RAFs. RETRACT the trade-off claim.\n", mean_adv);
else
printf(" => THE REVERSAL WAS A SIZE ARTIFACT. At matched realized size the two modes are\n"
" indistinguishable (%+.2f). Fission only looked better because propagule shrinks\n"
" compartments and smaller compartments hold smaller RAFs. There is no evidence here\n"
" that a bottleneck destroys closure — RETRACT the trade-off claim as stated.\n", mean_adv);
return 0;
}