import numpy as np import statistics as stats import multiprocessing as mp import csv, math, json, time from pathlib import Path # Independent replication of Rohde, Olson & Chang (Nature 431, 2004) # simple graph simulations (Tables 1 and 2). # # Operationalization of "exchange one pair of migrants per edge per generation": # for each undirected edge i--j and each generation, exactly one current # individual resident in i is treated as a migrant originating in j and exactly # one resident in j as a migrant originating in i. Both parents of a migrant are # sampled uniformly from its origin subpopulation; both parents of a non-migrant # are sampled uniformly from its resident subpopulation. Parent choices are with # replacement and independent. This is an independent implementation, not the # authors' original source code. GRAPHS = { 'One node': [[]], 'Three fully connected nodes': [[1,2],[0,2],[0,1]], 'Five fully connected nodes': [[j for j in range(5) if j != i] for i in range(5)], # Figure 1 edges: 1-2, 2-3, 3-4, 3-8, 4-5, 5-6, 5-7, 8-9, 8-10 (paper numbering) 'Ten-node graph shown in Fig. 1': [[1],[0,2],[1,3,7],[2,4],[3,5,6],[4],[4],[2,8,9],[7],[7]], } PUBLISHED_T = { 'One node': {1000:(10.8,0.4),2000:(11.8,0.4),4000:(12.8,0.4),8000:(13.9,0.3),16000:(14.8,0.4)}, 'Three fully connected nodes': {1000:(14.0,0.7),2000:(15.6,0.7),4000:(17.1,0.9),8000:(18.9,0.8),16000:(20.3,1.0)}, 'Five fully connected nodes': {1000:(14.0,0.5),2000:(15.8,0.5),4000:(17.8,0.5),8000:(19.6,0.5),16000:(21.5,0.6)}, 'Ten-node graph shown in Fig. 1': {1000:(21.1,1.3),2000:(24.3,1.5),4000:(27.6,1.5),8000:(30.5,1.5),16000:(33.8,1.7)}, } PUBLISHED_U = { 'One node': {1000:(20.8,1.6),2000:(22.6,1.5),4000:(24.6,1.5),8000:(26.5,1.6),16000:(28.3,1.4)}, 'Three fully connected nodes': {1000:(27.4,1.5),2000:(30.3,1.4),4000:(33.4,1.5),8000:(36.2,1.7),16000:(38.9,1.5)}, 'Five fully connected nodes': {1000:(25.9,1.3),2000:(28.9,1.4),4000:(32.1,1.7),8000:(35.3,1.5),16000:(37.9,1.4)}, 'Ten-node graph shown in Fig. 1': {1000:(46.3,2.7),2000:(53.0,2.7),4000:(59.8,2.7),8000:(66.8,2.9),16000:(73.6,2.7)}, } THEORY_INCREMENT_T = {'One node':1.00,'Three fully connected nodes':1.50,'Five fully connected nodes':1.75,'Ten-node graph shown in Fig. 1':3.00} THEORY_INCREMENT_U = {'One node':1.77,'Three fully connected nodes':2.77,'Five fully connected nodes':2.77,'Ten-node graph shown in Fig. 1':6.77} def split_sizes(N, G): q, r = divmod(N, G) return [q + (1 if i < r else 0) for i in range(G)] def run_once(args): graph_name, N, seed = args adj = GRAPHS[graph_name] G = len(adj) sizes = split_sizes(N, G) offsets = np.cumsum([0] + sizes) rng = np.random.default_rng(seed) # Bit k is set iff present-day individual k descends from this past individual. descendants = [1 << i for i in range(N)] everyone = (1 << N) - 1 T = None max_generations = 140 for generation in range(1, max_generations + 1): previous = [0] * N # For each node, mark exactly one migrant child from each adjacent node. origins = [] for node in range(G): m = sizes[node] origin = np.full(m, node, dtype=np.int16) if adj[node]: slots = np.arange(m) rng.shuffle(slots) for k, neigh in enumerate(adj[node]): origin[slots[k]] = neigh origins.append(origin) # Sample two independent parents for every child, both from its origin node. for node in range(G): lo, hi = int(offsets[node]), int(offsets[node+1]) m = sizes[node] origin = origins[node] p1 = np.empty(m, dtype=np.int32) p2 = np.empty(m, dtype=np.int32) for o in np.unique(origin): mask = (origin == o) count = int(mask.sum()) olo, ohi = int(offsets[o]), int(offsets[o+1]) p1[mask] = rng.integers(olo, ohi, size=count) p2[mask] = rng.integers(olo, ohi, size=count) for local_idx, child in enumerate(range(lo, hi)): bits = descendants[child] previous[int(p1[local_idx])] |= bits previous[int(p2[local_idx])] |= bits if T is None and any(bits == everyone for bits in previous): T = generation if T is not None and all(bits == 0 or bits == everyone for bits in previous): return graph_name, N, seed, T, generation descendants = previous raise RuntimeError(f'No IA point by {max_generations} generations: {graph_name}, N={N}, seed={seed}') def graph_metrics(adj): G = len(adj) if G == 1: return {'G':1,'R':0,'D':0,'Delta':None} dist=[] for s in range(G): d=[10**9]*G; d[s]=0; queue=[s] for x in queue: for y in adj[x]: if d[y] > d[x] + 1: d[y] = d[x] + 1; queue.append(y) dist.append(d) ecc=[max(row) for row in dist] R=min(ecc); D=max(ecc) centers=[i for i,e in enumerate(ecc) if e==R] # brute force minimal H_i among neighbors satisfying all nodes within R-1 of {i}+S import itertools Hs=[] for i in centers: neigh=adj[i] best=None for h in range(1,len(neigh)+1): for combo in itertools.combinations(neigh,h): A=(i,)+combo if all(min(dist[a][j] for a in A) <= R-1 for j in range(G)): best=h; break if best is not None: break Hs.append(best) H=min(Hs) Delta=(H-1)/H return {'G':G,'R':R,'D':D,'H':H,'Delta':Delta,'centers':[x+1 for x in centers]} def main(): reps=100 Ns=[1000,2000,4000,8000,16000] tasks=[] for gi, graph_name in enumerate(GRAPHS): for N in Ns: for r in range(reps): seed = 202608310000 + gi*1000000 + N*100 + r tasks.append((graph_name,N,seed)) t0=time.time() workers=min(12, max(1, mp.cpu_count()-1)) with mp.Pool(workers) as pool: raw=list(pool.imap_unordered(run_once,tasks,chunksize=1)) elapsed=time.time()-t0 grouped={} for graph_name,N,seed,T,U in raw: grouped.setdefault((graph_name,N),[]).append((T,U)) out=[] for graph_name in GRAPHS: for N in Ns: vals=grouped[(graph_name,N)] Ts=[x[0] for x in vals]; Us=[x[1] for x in vals] tm=stats.mean(Ts); tsd=stats.stdev(Ts) um=stats.mean(Us); usd=stats.stdev(Us) pt,ptsd=PUBLISHED_T[graph_name][N] pu,pusd=PUBLISHED_U[graph_name][N] out.append({ 'graph':graph_name,'N':N,'replicates':reps, 'T_rep_mean':tm,'T_rep_sd':tsd,'T_pub_mean':pt,'T_pub_sd':ptsd,'T_mean_diff':tm-pt, 'U_rep_mean':um,'U_rep_sd':usd,'U_pub_mean':pu,'U_pub_sd':pusd,'U_mean_diff':um-pu, }) csv_path=Path('/mnt/data/rohde_full_validation_results.csv') with csv_path.open('w',newline='') as f: w=csv.DictWriter(f,fieldnames=out[0].keys()); w.writeheader(); w.writerows(out) # raw per-run results for independent inspection raw_path=Path('/mnt/data/rohde_full_validation_raw_runs.csv') with raw_path.open('w',newline='') as f: w=csv.writer(f); w.writerow(['graph','N','seed','T','U']); w.writerows(sorted(raw,key=lambda x:(x[0],x[1],x[2]))) # Calculate doubling increments from replicate means and RMSE vs published tables. summary={} for graph_name in GRAPHS: rows=[r for r in out if r['graph']==graph_name] rows=sorted(rows,key=lambda r:r['N']) t_inc=[rows[i+1]['T_rep_mean']-rows[i]['T_rep_mean'] for i in range(4)] u_inc=[rows[i+1]['U_rep_mean']-rows[i]['U_rep_mean'] for i in range(4)] summary[graph_name]={ 'graph_metrics':graph_metrics(GRAPHS[graph_name]), 'T_theory_increment':THEORY_INCREMENT_T[graph_name], 'T_rep_increment_mean':stats.mean(t_inc), 'T_rep_increments':t_inc, 'U_theory_increment':THEORY_INCREMENT_U[graph_name], 'U_rep_increment_mean':stats.mean(u_inc), 'U_rep_increments':u_inc, 'T_RMSE_vs_published':math.sqrt(stats.mean([(r['T_mean_diff'])**2 for r in rows])), 'U_RMSE_vs_published':math.sqrt(stats.mean([(r['U_mean_diff'])**2 for r in rows])), } summary['_meta']={'replicates_per_cell':reps,'workers':workers,'elapsed_seconds':elapsed,'seed_scheme':'fixed deterministic seeds beginning 202608310000'} Path('/mnt/data/rohde_full_validation_summary.json').write_text(json.dumps(summary,indent=2)) print(f'Completed {len(tasks)} runs in {elapsed:.1f}s with {workers} workers') for r in out: print(f"{r['graph'][:10]:10s} N={r['N']:5d} T {r['T_rep_mean']:.2f} ({r['T_rep_sd']:.2f}) pub {r['T_pub_mean']:.1f}; U {r['U_rep_mean']:.2f} ({r['U_rep_sd']:.2f}) pub {r['U_pub_mean']:.1f}") print('\nDoubling checks:') for g,s in summary.items(): if g=='_meta': continue print(g, 'T avg incr', round(s['T_rep_increment_mean'],3), 'theory',s['T_theory_increment'], 'U avg incr',round(s['U_rep_increment_mean'],3),'theory',s['U_theory_increment']) if __name__=='__main__': main()