#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ founder-effect.py Science Journaling Club, Volume 1 Issue 1, Fall 2024, "Populations and Chance". THE QUESTION ------------ When a handful of individuals colonise somewhere new, they carry only a sample of the variation that existed in the population they left. We want to know three things about that sample. How much of the original variation does it hold? How much of what it does not hold is gone for good? And if somebody walked onto the island a few hundred generations later with a sequencer, could they still tell that the founding had happened? WHAT THIS PROGRAM IS -------------------- This is a computation, not an observation. The club has no field site, no island, no sampling permit, no sequencing machine and no organisms of any kind. Nothing in this file was measured in the physical world. Every number below comes out of pseudo-random sampling from a Wright-Fisher Markov chain that we wrote ourselves, seeded once at the top so that anybody who runs the file again gets exactly these digits. The Monte Carlo run is the experiment. Where the output says "measured" it means "estimated from our own simulated replicates", and nothing else. THE MODEL --------- One source population, L independent neutral loci, K_A = 20 possible allele states per locus (the K-allele mutation model of Crow & Kimura). The source is treated as an effectively infinite gamete pool with a fixed, known allele frequency spectrum p, generated once by GEM / Poisson-Dirichlet stick breaking at theta = 1.6 and printed in full in Section 2. Alleles below one copy in 2K are dropped and the rest renormalised, so every allele in the source is one that a population of size K could actually carry. FOUNDING. A founding group of k diploid individuals is formed by the random union of 2k gametes drawn from the source pool. Allele counts in the founding group are Multinomial(2k, p), independently at every locus. GROWTH. The new population grows from k toward carrying capacity K by capped geometric growth, N_{t+1} = min(K, round(N_t (1 + r))), with N held at a floor of 2. r is the per-generation growth rate and is swept. REPRODUCTION. Non-overlapping generations. Gametes mutate under the K-allele model at rate mu per copy per generation, so pool frequencies become q_i = p_i (1 - mu) + (mu / (K_A - 1)) (1 - p_i), and the next generation's counts are Multinomial(2 N_{t+1}, q). DIVERSITY. Expected heterozygosity, also called gene diversity, is computed as the plug-in quantity H = 1 - sum_i x_i^2 averaged over loci, where x is the allele frequency vector in the population. Allele count A is the number of states present at non-zero frequency. WHAT WE COMPARE AGAINST, AND WHY THESE ARE NOT APPROXIMATIONS ------------------------------------------------------------- 1. Founding retention. If X ~ Multinomial(2k, p), then E[sum_i (X_i/2k)^2] = sum_i p_i^2 + (1 / 2k) (1 - sum_i p_i^2), so E[H_new] = H_old (1 - 1/(2k)). This is exact for the model, not a large-k approximation. It is the closed form the study plan asks us to check, and Section 2 of the output checks it at fifteen founding sizes. 2. Allele loss. Allele i is absent from the founding group with probability exactly (1 - p_i)^{2k}, so the expected number of source alleles lost at founding is sum_i (1 - p_i)^{2k}. Section 4 checks this in aggregate and allele by allele. 3. Finite source. If the source is a real population of N diploids and the 2k founder gametes are drawn without replacement, the hypergeometric variance gives E[H_new] = H_old [ 1 - (1/(2k)) (2N - 2k)/(2N - 1) ]. At k = N the bracket is exactly 1, so nothing whatever is lost. Section 5 checks this, including the k = N case where the answer must be exact in every single replicate rather than on average. 4. The whole trajectory. Writing F = sum_i x_i^2 for homozygosity and a = 1 - mu K_A/(K_A - 1), b = 2 a mu/(K_A - 1) + K_A mu^2/(K_A - 1)^2, the mutation step gives E[sum q_i^2] = a^2 F + b exactly, and the multinomial step gives E[F_{t+1}] = E[sum q_i^2] (1 - 1/(2N_{t+1})) + 1/(2N_{t+1}). Chaining those two lines is an exact recursion for expected homozygosity at every generation of the recovery, for any size trajectory. Section 7 runs it beside the simulation. With mu = 0 it collapses to the familiar E[H_{t+1}] = E[H_t](1 - 1/(2N_{t+1})). 5. Detectability. Section 9 builds the club's own version of the heterozygosity excess test of Cornuet & Luikart (1996). After a bottleneck, rare alleles disappear faster than heterozygosity does, so the population carries more diversity than its own allele count would lead you to expect. We calibrate what "expect" means empirically, from control populations that were founded by the full carrying capacity and so lost essentially nothing, and we set the critical value from the controls so the false positive rate is 5 per cent by construction rather than by assumption. ASSUMPTIONS, STATED PLAINLY --------------------------- * Neutrality. No locus affects fitness. No selection of any kind, so nothing in this file speaks to adaptation on islands. * Loci are unlinked and independent. No recombination is modelled because none is needed; each locus is its own chain. * Random mating in one panmictic pool every generation. No assortative mating, no selfing beyond what random union of gametes implies, no spatial structure inside the new population. * Non-overlapping generations, everybody reproduces at once and then dies. No age structure, no overlapping cohorts, no seed bank or resting stage. * The census size is the effective size. Real populations have Ne below N, usually well below, so a real founder event of k individuals loses more than our k does. * Growth is deterministic given r. No demographic stochasticity in the size trajectory, no extinction of the new population, no environmental variance. The companion study on extinction by bad luck handles the case where the new population simply dies. * One-way colonisation. No migration back from the source and no further arrivals, so the only route back to diversity is mutation. * K-allele mutation with K_A = 20 states, symmetric, rate mu per copy per generation. Real microsatellites mutate stepwise and real sequence data has an effectively infinite allele space; both differ from this. * The source pool is infinite in Sections 3, 4, 6, 7 and 9, and finite only in Section 5 where finiteness is the point. LIMITATIONS ----------- Every one of those assumptions is wrong about some real colonisation and several are wrong about most. The strongest is the neutrality assumption: a real founding group on a real island meets a new selective environment immediately, and some of the variation it lost was variation it did not need. The second strongest is the deterministic growth: a founding group of two that grows at five per cent a generation is, in the real world, usually a founding group of two that dies. What the file does establish is what the sampling arithmetic alone implies, and that part is not a modelling choice. A group of k founders draws 2k gene copies, and 2k gene copies cannot hold more than 2k alleles no matter what happens afterwards. RUNTIME ------- Between four and eight minutes on a laptop, depending on what else it is doing. The run behind the published output took 276 seconds. Requires numpy. python analysis/founder-effect.py > analysis/founder-effect-output.txt """ from __future__ import annotations import math import sys import time import numpy as np # --------------------------------------------------------------------------- # The seed. Stated in the article. Change it and every number below moves a # little, inside the quoted Monte Carlo error. # --------------------------------------------------------------------------- MASTER_SEED = 20241108 # ---- the genetic model ----------------------------------------------------- K_ALLELES = 20 # allele states per locus, K-allele mutation model THETA_SPEC = 1.6 # stick-breaking parameter for the source spectrum K_CAP = 1000 # carrying capacity of the new population, diploids MU = THETA_SPEC / (4.0 * K_CAP) # = 4e-4 per copy per generation N_LOCI_TOTAL = 20 # loci in the source spectrum L_FOUND = 20 # loci used in the founding sweep L_DYN = 5 # loci used in the recovery grid L_DET = 15 # loci used in the detection experiment # ---- Section 3 and 4: the founding sweep ----------------------------------- FOUND_SIZES = [2, 3, 5, 8, 10, 12, 20, 30, 50, 75, 100, 150, 200, 350, 500] REPS_FOUND = 20_000 CHUNK_FOUND = 2_500 # ---- Section 5: finite source ---------------------------------------------- FINITE_N = 60 FINITE_K = [2, 5, 10, 20, 30, 45, 60] REPS_FINITE = 20_000 # ---- Section 7: the recovery grid ------------------------------------------ GRID_K = [2, 5, 10, 25, 50, 100, 250, 500] GRID_R = [0.05, 0.15, 0.40, 1.00] REPS_GRID = 5_000 GRID_TMIN = 50 GRID_TMAX = 170 GRID_SETTLE = 20 # generations held at carrying capacity before stopping # ---- Section 9: the detection window --------------------------------------- DET_K_ISL = 400 # island carrying capacity DET_R = 0.25 DET_FOUNDERS = [5, 25] DET_CONTROL = DET_K_ISL DET_REPS = 5_000 DET_T = 200 DET_SAMPLE = 60 # gene copies genotyped per locus, 30 diploids SENS_T = 120 # generations for the Section 8 sensitivity runs DET_CHECKS = [10, 13, 17, 22, 28, 36, 46, 60, 78, 100, 128, 160, 200] W = 118 # output rule width # =========================================================================== # small helpers # =========================================================================== def rule(ch="="): return ch * W def banner(title, ch="="): return "\n".join([rule(ch), title, rule(ch)]) def fmt(x, nd=6): return f"{x:.{nd}f}" def sigmas(obs, pred, se): """Standardised difference. A standard error of zero means every replicate gave the same answer, which is an exact result rather than a disagreement.""" if se <= 1e-12: return float("nan") return (obs - pred) / se def verdict(z, tol=3.0): if not np.isfinite(z): return "exact" return "agree" if abs(z) <= tol else "DISAGREE" # =========================================================================== # the source population # =========================================================================== def stick_break_spectrum(rng, theta, n_slots, min_freq): """ One locus of the source spectrum. GEM / Poisson-Dirichlet stick breaking at parameter theta, which is the neutral infinite-alleles equilibrium spectrum, truncated to n_slots states. Alleles rarer than min_freq are dropped and the remainder renormalised, because an allele below one copy in the population cannot be in the population. """ w = np.zeros(n_slots, dtype=np.float64) rem = 1.0 for i in range(n_slots - 1): v = rng.beta(1.0, theta) w[i] = rem * v rem *= (1.0 - v) w[n_slots - 1] = rem w = np.sort(w)[::-1] w[w < min_freq] = 0.0 tot = w.sum() if tot <= 0: w[0] = 1.0 tot = 1.0 return w / tot def build_source(rng): min_freq = 1.0 / (2.0 * K_CAP) spec = np.zeros((N_LOCI_TOTAL, K_ALLELES), dtype=np.float64) for l in range(N_LOCI_TOTAL): spec[l] = stick_break_spectrum(rng, THETA_SPEC, K_ALLELES, min_freq) return spec # =========================================================================== # population dynamics # =========================================================================== def size_path_geometric(k, r, K, tmin, tmax, settle): """Capped geometric growth. Returns N_1 .. N_T, the sizes AFTER founding.""" sizes = [] N = float(k) at_cap = 0 t = 0 while True: N = min(float(K), N * (1.0 + r)) n = int(round(N)) n = max(2, min(K, n)) sizes.append(n) t += 1 if n >= K: at_cap += 1 if (t >= tmin and at_cap >= settle) or t >= tmax: break return np.array(sizes, dtype=np.int64) def size_path_logistic(k, r, K, tmin, tmax, settle): """Discrete logistic growth, used only for the Section 8 sensitivity check.""" sizes = [] N = float(k) at_cap = 0 t = 0 while True: N = min(float(K), max(1.0, N + r * N * (1.0 - N / K))) n = int(round(N)) n = max(2, min(K, n)) sizes.append(n) t += 1 if n >= K: at_cap += 1 if (t >= tmin and at_cap >= settle) or t >= tmax: break return np.array(sizes, dtype=np.int64) def mutate(p, mu, ka): """K-allele mutation applied to a frequency array of shape (rows, ka).""" if mu <= 0.0: return p.copy() return p * (1.0 - mu) + (mu / (ka - 1.0)) * (1.0 - p) def recursion_F(F0, sizes, mu, ka): """ Exact recursion for expected homozygosity along a size path. Returns array of length len(sizes)+1, element 0 being F0. """ a = 1.0 - mu * ka / (ka - 1.0) b = 2.0 * a * mu / (ka - 1.0) + ka * mu * mu / ((ka - 1.0) ** 2) out = np.empty(len(sizes) + 1, dtype=np.float64) out[0] = F0 F = F0 for i, n in enumerate(sizes): G = a * a * F + b F = G * (1.0 - 1.0 / (2.0 * n)) + 1.0 / (2.0 * n) out[i + 1] = F return out def evolve(rng, counts, sizes, mu, ka, record=None): """ Run counts (rows, ka) forward along the size path. counts rows are (replicate, locus) pairs already flattened. If record is a dict mapping generation index -> callback, the callback is handed the counts array. Generation 0 is the founding group itself. """ if record is not None and 0 in record: record[0](counts) for t, n in enumerate(sizes, start=1): tot = counts.sum(axis=1, keepdims=True).astype(np.float64) p = counts / tot q = mutate(p, mu, ka) # shrink by one part in 1e12 so floating point can never let the leading # categories sum above 1, which numpy's multinomial refuses q /= (q.sum(axis=1, keepdims=True) * (1.0 + 1e-12)) counts = rng.multinomial(2 * int(n), q).astype(np.int32) if record is not None and t in record: record[t](counts) return counts def het_and_alleles(counts, n_reps, n_loci): """ counts has shape (n_reps * n_loci, ka). Returns per-replicate mean gene diversity across loci, and per-replicate mean allele count across loci. """ tot = counts.sum(axis=1).astype(np.float64) f = counts.astype(np.float64) / tot[:, None] h = 1.0 - (f * f).sum(axis=1) a = (counts > 0).sum(axis=1).astype(np.float64) return (h.reshape(n_reps, n_loci).mean(axis=1), a.reshape(n_reps, n_loci).mean(axis=1)) # =========================================================================== # main # =========================================================================== def main(): t_start = time.time() ss = np.random.SeedSequence(MASTER_SEED) streams = ss.spawn(12) rng_spec = np.random.default_rng(streams[0]) out = [] P = out.append P(rule()) P("EVERYTHING A SMALL FOUNDING GROUP LOSES ON THE WAY") P("Monte Carlo study of the founder effect: sampling, recovery and detectability") P("Science Journaling Club, Volume 1 Issue 1, Fall 2024") P(rule()) P("") P("This output is generated entirely by simulation. No organism, population,") P("island, field site or laboratory measurement is involved anywhere in it.") P("") P(f"numpy version : {np.__version__}") P(f"python version : {sys.version.split()[0]}") P(f"master seed : {MASTER_SEED}") P("generator : numpy PCG64, one spawned stream per experiment") P(f"allele states per locus : {K_ALLELES} (K-allele mutation model)") P(f"mutation rate mu : {MU:.3e} per copy per generation") P(f"carrying capacity K : {K_CAP:,} diploids") P(f"theta = 4 K mu : {4.0 * K_CAP * MU:.4f}") P(f"loci in source spectrum : {N_LOCI_TOTAL}") P("") P("Founding: 2k gametes drawn from the source pool, Multinomial(2k, p).") P("Growth: N_{t+1} = min(K, round(N_t (1 + r))).") P("Breeding: gametes mutate under the K-allele model, then Multinomial(2N, q).") # ----------------------------------------------------------------------- # Section 1, the source # ----------------------------------------------------------------------- spec = build_source(rng_spec) src_hom = (spec * spec).sum(axis=1) src_het = 1.0 - src_hom src_alleles = (spec > 0).sum(axis=1) H_SRC = float(src_het.mean()) H_SRC_FOUND = float(src_het[:L_FOUND].mean()) H_SRC_DYN = float(src_het[:L_DYN].mean()) H_SRC_DET = float(src_het[:L_DET].mean()) A_SRC = float(src_alleles.mean()) P("") P(banner("SECTION 1. THE SOURCE POPULATION")) P("") P("The spectrum below is drawn once, from stream 0, and then held fixed for the") P("rest of the study. It is the known quantity the founding groups sample from.") P("") P(" locus alleles gene diversity commonest rarest") P(" " + "-" * 56) for l in range(N_LOCI_TOTAL): nz = spec[l][spec[l] > 0] P(f" {l:5d} {src_alleles[l]:7d} {src_het[l]:14.6f} " f"{nz.max():9.5f} {nz.min():8.6f}") P(" " + "-" * 56) P(f" mean {A_SRC:7.2f} {H_SRC:14.6f}") P("") P(f"H_source over all {N_LOCI_TOTAL} loci : {H_SRC:.6f}") P(f"H_source over the {L_FOUND} founding-sweep loci : {H_SRC_FOUND:.6f}") P(f"H_source over the {L_DYN} recovery-grid loci : {H_SRC_DYN:.6f}") P(f"H_source over the {L_DET} detection loci : {H_SRC_DET:.6f}") P(f"mean alleles per locus in the source : {A_SRC:.3f}") P(f"total distinct alleles across all loci : {int(src_alleles.sum())}") # frequency classes of the source, for the spectrum figure later edges = np.array([0.0, 0.01, 0.02, 0.05, 0.10, 0.20, 0.40, 1.01]) flat_src = spec[:L_FOUND].ravel() flat_src = flat_src[flat_src > 0] src_hist = np.histogram(flat_src, bins=edges)[0].astype(np.float64) / L_FOUND P("") P("Allele frequency classes in the source, mean count per locus") P(" class alleles/locus") for i in range(len(edges) - 1): lo, hi = edges[i], min(edges[i + 1], 1.0) P(f" {lo:5.2f} to {hi:5.2f} {src_hist[i]:8.3f}") # ----------------------------------------------------------------------- # Section 2, the founding sweep # ----------------------------------------------------------------------- rng_found = np.random.default_rng(streams[1]) p_found = spec[:L_FOUND] p_tile = np.repeat(p_found[None, :, :], CHUNK_FOUND, axis=0).reshape(-1, K_ALLELES) src_present = (p_found > 0) src_alleles_found = int(src_present.sum()) P("") P(banner("SECTION 2. THE FOUNDING DRAW, AND THE CLOSED FORM")) P("") P("Prediction, exact for this model and not an approximation:") P(" E[H_new] = H_old (1 - 1/(2k)).") P(f"H_old here is {H_SRC_FOUND:.6f}, the mean gene diversity of the {L_FOUND} source loci.") P(f"Replicates per founding size: {REPS_FOUND:,}.") P("") hdr = (" k retention predicted diff SE z " "verdict H_new sim H_new pred") P(hdr) P(" " + "-" * (len(hdr) - 2)) found_rows = [] conv_store = {} for k in FOUND_SIZES: hs = np.empty(REPS_FOUND, dtype=np.float64) as_ = np.empty(REPS_FOUND, dtype=np.float64) lost = np.empty(REPS_FOUND, dtype=np.float64) per_allele_present = np.zeros((L_FOUND, K_ALLELES), dtype=np.int64) done = 0 while done < REPS_FOUND: m = min(CHUNK_FOUND, REPS_FOUND - done) pt = p_tile if m == CHUNK_FOUND else np.repeat( p_found[None, :, :], m, axis=0).reshape(-1, K_ALLELES) cnt = rng_found.multinomial(2 * k, pt).astype(np.int32) h, a = het_and_alleles(cnt, m, L_FOUND) hs[done:done + m] = h as_[done:done + m] = a pres = (cnt > 0).reshape(m, L_FOUND, K_ALLELES) per_allele_present += pres.sum(axis=0) lost[done:done + m] = (src_present[None, :, :] & ~pres).sum(axis=(1, 2)) done += m ret = hs / H_SRC_FOUND ret_mean = float(ret.mean()) ret_se = float(ret.std(ddof=1) / math.sqrt(REPS_FOUND)) pred = 1.0 - 1.0 / (2.0 * k) z = sigmas(ret_mean, pred, ret_se) P(f" {k:6d} {ret_mean:9.6f} {pred:10.6f} {ret_mean - pred:+9.6f} " f"{ret_se:7.6f} {z:+6.2f} {verdict(z):9s} {hs.mean():9.6f} " f"{pred * H_SRC_FOUND:10.6f}") exp_lost = float((np.power(1.0 - p_found, 2 * k) * src_present).sum()) found_rows.append(dict( k=k, ret=ret_mean, ret_se=ret_se, pred=pred, z=z, h=float(hs.mean()), h_se=float(hs.std(ddof=1) / math.sqrt(REPS_FOUND)), alleles=float(as_.mean()), alleles_se=float(as_.std(ddof=1) / math.sqrt(REPS_FOUND)), lost=float(lost.mean()), lost_se=float(lost.std(ddof=1) / math.sqrt(REPS_FOUND)), exp_lost=exp_lost, per_allele=per_allele_present.copy())) if k in (2, 10, 50): conv_store[k] = ret.copy() zs = np.array([r["z"] for r in found_rows]) P(" " + "-" * (len(hdr) - 2)) P(f" largest |z| over the {len(FOUND_SIZES)} founding sizes: {np.abs(zs).max():.2f}") P(f" mean z: {zs.mean():+.3f} (expected near 0, SD near 1; realised SD " f"{zs.std(ddof=1):.3f})") P("") P(" Every founding size agrees with H_old(1 - 1/(2k)) inside three standard") P(" errors. There is no free parameter in that comparison.") # ----------------------------------------------------------------------- # Section 3, alleles lost # ----------------------------------------------------------------------- P("") P(banner("SECTION 3. ALLELES LOST AT THE MOMENT OF FOUNDING")) P("") P("Allele i survives the draw unless all 2k gametes miss it, so") P(" P(lost) = (1 - p_i)^{2k}, E[alleles lost] = sum_i (1 - p_i)^{2k}.") P(f"The {L_FOUND} source loci carry {src_alleles_found} distinct alleles between them.") P("") hdr2 = (" k lost sim SE lost pred diff z verdict " "alleles/locus % of source") P(hdr2) P(" " + "-" * (len(hdr2) - 2)) for r in found_rows: z = sigmas(r["lost"], r["exp_lost"], r["lost_se"]) P(f" {r['k']:6d} {r['lost']:8.4f} {r['lost_se']:7.4f} " f"{r['exp_lost']:9.4f} {r['lost'] - r['exp_lost']:+8.4f} {z:+5.2f} " f"{verdict(z):7s} {r['alleles']:13.4f} {100.0 * r['alleles'] / A_SRC:10.2f}") P(" " + "-" * (len(hdr2) - 2)) zl = np.array([sigmas(r["lost"], r["exp_lost"], r["lost_se"]) for r in found_rows]) P(f" largest |z| on allele loss: {np.abs(zl).max():.2f}") P("") P("Allele by allele, at k = 2 and k = 50. Only the rarest and commonest few are") P("printed; the full set behaves the same way. p is the source frequency,") P("P(lost) the binomial prediction, and the last column the simulated rate.") P("") for kk in (2, 50): row = [r for r in found_rows if r["k"] == kk][0] pa = row["per_allele"] flat = [] for l in range(L_FOUND): for i in range(K_ALLELES): if spec[l, i] > 0: pl = 1.0 - pa[l, i] / REPS_FOUND flat.append((spec[l, i], (1.0 - spec[l, i]) ** (2 * kk), pl)) flat.sort() P(f" k = {kk}") P(" p P(lost) pred P(lost) sim diff z") P(" " + "-" * 58) show = flat[:6] + flat[len(flat) // 2 - 2:len(flat) // 2 + 2] + flat[-4:] for p_i, pr, ob in show: se = math.sqrt(max(pr * (1 - pr), 1e-12) / REPS_FOUND) P(f" {p_i:9.6f} {pr:12.6f} {ob:12.6f} {ob - pr:+9.6f} " f"{(ob - pr) / se:+5.2f}") allp = np.array([f[1] for f in flat]) allo = np.array([f[2] for f in flat]) P(f" across all {len(flat)} alleles: mean predicted {allp.mean():.6f}, " f"mean simulated {allo.mean():.6f}, mean |diff| {np.abs(allp - allo).mean():.6f}") P("") # ----------------------------------------------------------------------- # Section 4, finite source # ----------------------------------------------------------------------- P(banner("SECTION 4. A FINITE SOURCE, AND THE CASE k = N")) P("") P(f"The source is now a real population of N = {FINITE_N} diploids, {2 * FINITE_N} gene copies,") P("whose composition is drawn once from the spectrum of Section 1. Founder") P("gametes are taken without replacement, so") P(" E[H_new] = H_old [1 - (1/(2k)) (2N - 2k)/(2N - 1)].") P("At k = N the bracket is exactly 1 and the founding group is the source.") P("") rng_fin = np.random.default_rng(streams[2]) fin_counts = np.zeros((L_FOUND, K_ALLELES), dtype=np.int64) for l in range(L_FOUND): fin_counts[l] = rng_fin.multinomial(2 * FINITE_N, spec[l]) fin_f = fin_counts / (2.0 * FINITE_N) fin_h = 1.0 - (fin_f * fin_f).sum(axis=1) H_FIN = float(fin_h.mean()) A_FIN = int((fin_counts > 0).sum()) P(f"realised finite source: H_old = {H_FIN:.6f}, " f"{A_FIN} distinct alleles over {L_FOUND} loci") P("") hdr3 = (" k H sim SE H pred infinite H pred finite diff" " z verdict alleles lost") P(hdr3) P(" " + "-" * (len(hdr3) - 2)) finite_rows = [] for k in FINITE_K: h_rep = np.zeros((REPS_FINITE, L_FOUND), dtype=np.float64) lost_rep = np.zeros(REPS_FINITE, dtype=np.float64) for l in range(L_FOUND): draw = rng_fin.multivariate_hypergeometric( fin_counts[l], 2 * k, size=REPS_FINITE) f = draw / (2.0 * k) h_rep[:, l] = 1.0 - (f * f).sum(axis=1) lost_rep += ((fin_counts[l] > 0)[None, :] & (draw == 0)).sum(axis=1) hm = h_rep.mean(axis=1) h_mean = float(hm.mean()) h_se = float(hm.std(ddof=1) / math.sqrt(REPS_FINITE)) pred_inf = H_FIN * (1.0 - 1.0 / (2.0 * k)) corr = (2.0 * FINITE_N - 2.0 * k) / (2.0 * FINITE_N - 1.0) pred_fin = H_FIN * (1.0 - corr / (2.0 * k)) z = sigmas(h_mean, pred_fin, h_se) P(f" {k:6d} {h_mean:8.6f} {h_se:7.6f} {pred_inf:13.6f} " f"{pred_fin:13.6f} {h_mean - pred_fin:+8.6f} {z:+5.2f} " f"{verdict(z):7s} {lost_rep.mean():12.4f}") finite_rows.append((k, h_mean, h_se, pred_inf, pred_fin, z, float(lost_rep.mean()), float(lost_rep.max()))) P(" " + "-" * (len(hdr3) - 2)) kN = finite_rows[-1] P("") P(f" At k = N = {FINITE_N}: simulated H = {kN[1]:.10f}, source H = {H_FIN:.10f},") P(f" difference {kN[1] - H_FIN:+.2e}. Alleles lost: mean {kN[6]:.4f}, " f"maximum over {REPS_FINITE:,} replicates {kN[7]:.0f}.") P(" Nothing is lost, in any replicate, exactly as the hypergeometric variance says.") P("") P(" Note the fourth column. Using the infinite-source formula on a finite source") P(f" overstates the loss badly once k approaches N: at k = 30 it predicts " f"{finite_rows[4][3]:.6f}") P(f" where the right answer is {finite_rows[4][4]:.6f}. The club got this wrong once.") # ----------------------------------------------------------------------- # Section 5, convergence # ----------------------------------------------------------------------- P("") P(banner("SECTION 5. CONVERGENCE OF THE MONTE CARLO ESTIMATE")) P("") P("Running mean of the retention ratio H_new / H_old against replicate count,") P("for three founding sizes, with the analytic value it is converging on.") P("") conv_points = [50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000] P(" reps k=2 running k=2 target k=10 running k=10 target" " k=50 running k=50 target") P(" " + "-" * 104) conv_table = {} for k in (2, 10, 50): r = conv_store[k] conv_table[k] = np.cumsum(r) / np.arange(1, len(r) + 1) for n in conv_points: if n > REPS_FOUND: continue line = f" {n:8d} " for k in (2, 10, 50): tgt = 1.0 - 1.0 / (2.0 * k) line += f" {conv_table[k][n - 1]:11.6f} {tgt:10.6f}" P(line) P(" " + "-" * 104) for k in (2, 10, 50): r = conv_store[k] tgt = 1.0 - 1.0 / (2.0 * k) se = r.std(ddof=1) / math.sqrt(len(r)) P(f" k = {k:3d}: final {r.mean():.6f}, target {tgt:.6f}, SE {se:.6f}, " f"|final - target| / SE = {abs(r.mean() - tgt) / se:.2f}") # ----------------------------------------------------------------------- # Section 6, the recovery grid # ----------------------------------------------------------------------- P("") P(banner("SECTION 6. THE RECOVERY GRID")) P("") P(f"{len(GRID_K)} founding sizes x {len(GRID_R)} growth rates = " f"{len(GRID_K) * len(GRID_R)} cells, {REPS_GRID:,} replicates each,") P(f"{L_DYN} loci per replicate, carrying capacity {K_CAP:,}.") P("Each cell runs until the population has sat at carrying capacity for") P(f"{GRID_SETTLE} generations, with a floor of {GRID_TMIN} and a ceiling of {GRID_TMAX}.") P("") P("H_end is gene diversity at the last simulated generation. The 'recursion'") P("column is the exact expectation from the two-line recursion in the docstring,") P("run along the same size path. No fitting of any kind.") P("") hdr4 = (" k r gens N_end H_found H_end recursion diff" " z verdict A_found A_end % H kept % A kept") P(hdr4) P(" " + "-" * (len(hdr4) - 2)) rng_grid = np.random.default_rng(streams[3]) p_dyn = spec[:L_DYN] p_dyn_tile = np.repeat(p_dyn[None, :, :], REPS_GRID, axis=0).reshape(-1, K_ALLELES) F_SRC_DYN = float((p_dyn * p_dyn).sum(axis=1).mean()) grid_rows = [] traj_store = {} for k in GRID_K: for r in GRID_R: sizes = size_path_geometric(k, r, K_CAP, GRID_TMIN, GRID_TMAX, GRID_SETTLE) cnt = rng_grid.multinomial(2 * k, p_dyn_tile).astype(np.int32) h0, a0 = het_and_alleles(cnt, REPS_GRID, L_DYN) keep = k in (2, 25, 250) traj = {"h": [float(h0.mean())], "a": [float(a0.mean())], "n": [k], "t": [0]} if keep else None rec = {} if keep: def make_cb(store): def cb(c): hh, aa = het_and_alleles(c, REPS_GRID, L_DYN) store["h"].append(float(hh.mean())) store["a"].append(float(aa.mean())) return cb cb = make_cb(traj) for t in range(1, len(sizes) + 1): rec[t] = cb cnt = evolve(rng_grid, cnt, sizes, MU, K_ALLELES, record=rec if keep else None) h1, a1 = het_and_alleles(cnt, REPS_GRID, L_DYN) h_end = float(h1.mean()) h_se = float(h1.std(ddof=1) / math.sqrt(REPS_GRID)) F0 = F_SRC_DYN + (1.0 - F_SRC_DYN) / (2.0 * k) Fpath = recursion_F(F0, sizes, MU, K_ALLELES) h_pred = 1.0 - Fpath[-1] z = sigmas(h_end, h_pred, h_se) if keep: traj["n"] = [k] + list(sizes) traj["t"] = list(range(0, len(sizes) + 1)) traj["pred"] = list(1.0 - Fpath) traj_store[(k, r)] = traj P(f" {k:6d} {r:5.2f} {len(sizes):5d} {sizes[-1]:5d} " f"{h0.mean():9.6f} {h_end:8.6f} {h_pred:10.6f} " f"{h_end - h_pred:+8.6f} {z:+5.2f} {verdict(z):7s} " f"{a0.mean():7.3f} {a1.mean():6.3f} " f"{100.0 * h_end / H_SRC_DYN:7.2f} {100.0 * a1.mean() / A_SRC:7.2f}") grid_rows.append(dict( k=k, r=r, gens=len(sizes), h0=float(h0.mean()), h0_se=float(h0.std(ddof=1) / math.sqrt(REPS_GRID)), h_end=h_end, h_se=h_se, h_pred=h_pred, z=z, a0=float(a0.mean()), a_end=float(a1.mean()), a_end_se=float(a1.std(ddof=1) / math.sqrt(REPS_GRID)))) sys.stderr.write(f" grid k={k} done at {time.time() - t_start:.0f}s\n") sys.stderr.flush() P(" " + "-" * (len(hdr4) - 2)) zg = np.array([g["z"] for g in grid_rows]) P(f" largest |z| against the exact recursion over " f"{len(grid_rows)} cells: {np.abs(zg).max():.2f}") P(f" mean z {zg.mean():+.3f}, SD {zg.std(ddof=1):.3f}") # ----------------------------------------------------------------------- # Section 7, the long view # ----------------------------------------------------------------------- P("") P(banner("SECTION 7. HOW LONG UNTIL IT COMES BACK")) P("") P("The recursion is exact and costs nothing to run, so once Section 6 has shown") P("it tracks the simulation we can push it far past what we can afford to") P("simulate. Below, a population founded by k individuals grows at r = 0.40 to") P(f"K = {K_CAP:,} and then sits there. Gene diversity as a percentage of the source,") P("at generations we could never reach by Monte Carlo.") P("") long_gens = [0, 1, 5, 10, 25, 50, 100, 250, 500, 1000, 2000, 5000, 10000, 20000] P(" k" + "".join(f"{g:>9d}" for g in long_gens)) P(" " + "-" * (8 + 9 * len(long_gens))) long_rows = {} for k in [2, 5, 10, 25, 50, 100, 250, 500]: sizes = size_path_geometric(k, 0.40, K_CAP, 1, max(long_gens), 10 ** 9) F0 = F_SRC_DYN + (1.0 - F_SRC_DYN) / (2.0 * k) Fp = recursion_F(F0, sizes, MU, K_ALLELES) hp = 1.0 - Fp vals = [100.0 * hp[min(g, len(hp) - 1)] / H_SRC_DYN for g in long_gens] long_rows[k] = vals P(f" {k:6d}" + "".join(f"{v:9.2f}" for v in vals)) P(" " + "-" * (8 + 9 * len(long_gens))) # equilibrium of the recursion at N = K a = 1.0 - MU * K_ALLELES / (K_ALLELES - 1.0) b = 2.0 * a * MU / (K_ALLELES - 1.0) + K_ALLELES * MU * MU / ((K_ALLELES - 1.0) ** 2) c1 = (1.0 - 1.0 / (2.0 * K_CAP)) F_eq = (b * c1 + 1.0 / (2.0 * K_CAP)) / (1.0 - a * a * c1) H_eq = 1.0 - F_eq P(f" mutation-drift equilibrium of this model at N = K: H = {H_eq:.6f}, " f"{100.0 * H_eq / H_SRC_DYN:.2f}% of the source") P(" The source spectrum itself sits slightly above that equilibrium, which is why") P(" even the unbottlenecked ceiling in the table is below 100 per cent.") P("") half = {} P(" Approach to the model's own equilibrium. The last column is the first") P(" generation from which expected gene diversity stays inside one per cent of") P(" H_eq for good, coming from whichever side it started on.") P("") P(" k H at founding minimum H at gen half the climb back within 1%") P(" " + "-" * 82) for k in [2, 5, 10, 25, 50, 100, 250, 500]: sizes = size_path_geometric(k, 0.40, K_CAP, 1, 60000, 10 ** 9) F0 = F_SRC_DYN + (1.0 - F_SRC_DYN) / (2.0 * k) hp = 1.0 - recursion_F(F0, sizes, MU, K_ALLELES) imin = int(np.argmin(hp)) if hp[imin] < H_eq - 1e-9: target = H_eq - 0.5 * (H_eq - hp[imin]) idx = np.where(hp[imin:] >= target)[0] hg = int(idx[0]) + imin if len(idx) else -1 else: hg = -1 outside = np.where(np.abs(hp - H_eq) > 0.01 * H_eq)[0] ng = int(outside[-1]) + 1 if len(outside) else 0 half[k] = hg dips = hp[imin] < H_eq - 1e-9 P(f" {k:6d} {hp[0]:14.6f} {hp[imin]:11.6f} " f"{(str(imin) if dips else '-'):>6s} " f"{(str(hg) if hg >= 0 else 'never dips below'):>19s} {ng:9d}") P(" " + "-" * 82) # ----------------------------------------------------------------------- # Section 8, sensitivity to two modelling choices # ----------------------------------------------------------------------- P("") P(banner("SECTION 8. WHERE A DIFFERENT MODELLING CHOICE CHANGES THE ANSWER")) P("") P("Two choices we made are arbitrary. Both are re-run here for k = 5, r = 0.15,") P(f"{REPS_GRID:,} replicates, and all three variants are run for exactly") P(f"{SENS_T} generations so that only the choice under test differs.") P("") rng_sens = np.random.default_rng(streams[4]) sens = {} for label, path_fn, mu in [ ("capped geometric growth, mutation on", size_path_geometric, MU), ("discrete logistic growth, mutation on", size_path_logistic, MU), ("capped geometric growth, mutation off", size_path_geometric, 0.0)]: sizes = path_fn(5, 0.15, K_CAP, SENS_T, SENS_T, 10 ** 9) cnt = rng_sens.multinomial(2 * 5, p_dyn_tile).astype(np.int32) cnt = evolve(rng_sens, cnt, sizes, mu, K_ALLELES) h1, a1 = het_and_alleles(cnt, REPS_GRID, L_DYN) F0 = F_SRC_DYN + (1.0 - F_SRC_DYN) / (2.0 * 5) hp = 1.0 - recursion_F(F0, sizes, mu, K_ALLELES)[-1] se = float(h1.std(ddof=1) / math.sqrt(REPS_GRID)) sens[label] = (len(sizes), float(h1.mean()), se, hp, float(a1.mean()), int(sizes.sum())) P(f" {label:40s} gens {len(sizes):4d} H_end {h1.mean():.6f} " f"+/- {se:.6f} recursion {hp:.6f} A_end {a1.mean():.3f}") g = sens["capped geometric growth, mutation on"] lg = sens["discrete logistic growth, mutation on"] nm = sens["capped geometric growth, mutation off"] P("") P(f" Individual-generations lived over the {SENS_T} generations: {g[5]:,} under capped") P(f" geometric growth against {lg[5]:,} under discrete logistic growth. Logistic") P(" growth holds the population below capacity for longer, and its final gene") P(f" diversity differs from the geometric case by {100.0 * (lg[1] - g[1]) / g[1]:+.2f}%, " f"which is {lg[1] - g[1]:+.6f}") P(f" against a standard error of {math.hypot(g[2], lg[2]):.6f} on the difference.") P(f" Switching mutation off changes final gene diversity by " f"{100.0 * (nm[1] - g[1]) / g[1]:+.2f}% and the") P(f" final allele count by {100.0 * (nm[4] - g[4]) / g[4]:+.2f}%, over the same " f"{SENS_T} generations.") # ----------------------------------------------------------------------- # Section 9, the detection window # ----------------------------------------------------------------------- P("") P(banner("SECTION 9. HOW LONG THE SIGNATURE STAYS VISIBLE")) P("") P("A new population on an island of carrying capacity " f"{DET_K_ISL}, growth rate {DET_R}.") P(f"{DET_REPS:,} replicates per cell, {L_DET} loci, and at each checkpoint a") P(f"genetic sample of {DET_SAMPLE} gene copies ({DET_SAMPLE // 2} diploids) per locus,") P("drawn with replacement from the population.") P("") P("Checkpoints start at generation 10 because before that the smaller founding") P(f"groups hold fewer than {DET_SAMPLE // 2} individuals, and a {DET_SAMPLE // 2}-individual sample is") P("not something anybody could take.") P("") P("The control is a population founded by the full island capacity. It lost") P("essentially nothing at the founding, but it lives in the same small place and") P("drifts at the same rate afterwards. The test therefore measures the extra") P("signature left by a SMALL founding group, over and above the signature left by") P("living somewhere small. That is a harder test than comparing against the source.") P("") P("The statistic is the heterozygosity excess of Cornuet & Luikart (1996). A") P("bottleneck kills rare alleles faster than it kills heterozygosity, so the") P("population carries more diversity than its own allele count implies. We") P("calibrate 'implies' from the controls, which were founded by the full") P(f"{DET_CONTROL} individuals and so lost essentially nothing, then count how many") P(f"of the {L_DET} loci sit above that control curve. The critical value is the") P("smallest count whose control tail probability is at or below 0.05, so the") P("false positive rate is 5 per cent by construction.") P("") rng_det = np.random.default_rng(streams[5]) p_det = spec[:L_DET] p_det_tile = np.repeat(p_det[None, :, :], DET_REPS, axis=0).reshape(-1, K_ALLELES) det_sizes = {} det_results = {} for k in DET_FOUNDERS + [DET_CONTROL]: sizes = size_path_geometric(k, DET_R, DET_K_ISL, DET_T, DET_T, 10 ** 9) det_sizes[k] = sizes cnt = rng_det.multinomial(2 * k, p_det_tile).astype(np.int32) store = {} def make_cb(store, gen): def cb(c): tot = c.sum(axis=1).astype(np.float64) f = c / tot[:, None] samp = rng_det.multinomial(DET_SAMPLE, f).astype(np.int32) x = samp / float(DET_SAMPLE) hh = (1.0 - (x * x).sum(axis=1)) * DET_SAMPLE / (DET_SAMPLE - 1.0) aa = (samp > 0).sum(axis=1) store[gen] = (hh.reshape(DET_REPS, L_DET), aa.reshape(DET_REPS, L_DET)) return cb rec = {} for gen in DET_CHECKS: if gen <= len(sizes): rec[gen] = make_cb(store, gen) evolve(rng_det, cnt, sizes, MU, K_ALLELES, record=rec) det_results[k] = store sys.stderr.write(f" detection k={k} done at {time.time() - t_start:.0f}s\n") sys.stderr.flush() # calibration and power P(" Control null: distribution of the number of loci above the control curve.") P("") hdr5 = (" gen N k=5 crit false pos power k=5 " "power k=25 H k=5 H ctrl A k=5 A ctrl") P(hdr5) P(" " + "-" * (len(hdr5) - 2)) det_power = {k: [] for k in DET_FOUNDERS} det_gens = [] det_extra = [] for gen in DET_CHECKS: if gen not in det_results[DET_CONTROL]: continue hc, ac = det_results[DET_CONTROL][gen] # control curve: mean sampled H at each observed allele count a_flat = ac.ravel() h_flat = hc.ravel() amax = int(a_flat.max()) curve = np.full(amax + 2, np.nan) for av in range(1, amax + 1): m = a_flat == av if m.sum() >= 30: curve[av] = h_flat[m].mean() idx = np.where(np.isfinite(curve))[0] if len(idx) >= 2: curve = np.interp(np.arange(len(curve)), idx, curve[idx]) else: curve = np.full(amax + 2, float(h_flat.mean())) def sign_count(h, a): exp = curve[np.clip(a, 0, len(curve) - 1)] return (h > exp).sum(axis=1) s_ctrl = sign_count(hc, ac) # critical value: smallest c with P(S >= c | control) <= 0.05 crit = L_DET + 1 for c in range(L_DET + 1): if (s_ctrl >= c).mean() <= 0.05: crit = c break fp = float((s_ctrl >= crit).mean()) n_at = int(det_sizes[5][gen - 1]) if gen > 0 else 5 row = [gen, n_at, crit, fp] for k in DET_FOUNDERS: hk, ak = det_results[k][gen] pw = float((sign_count(hk, ak) >= crit).mean()) det_power[k].append(pw) row.append(pw) h5 = det_results[5][gen][0].mean() a5 = det_results[5][gen][1].mean() P(f" {gen:6d} {row[1]:6d} {crit:5d} {fp:9.4f} " f"{row[4]:12.4f} {row[5]:10.4f} {h5:7.4f} {hc.mean():7.4f} " f"{a5:7.3f} {ac.mean():7.3f}") det_gens.append(gen) det_extra.append((gen, float(h5), float(hc.mean()), float(a5), float(ac.mean()))) P(" " + "-" * (len(hdr5) - 2)) P("") for k in DET_FOUNDERS: pw = np.array(det_power[k]) gg = np.array(det_gens) above50 = gg[pw >= 0.5] above80 = gg[pw >= 0.8] P(f" k = {k}: power stays at or above 0.80 through generation " f"{above80.max() if len(above80) else 'never reached'}, and at or above") P(f" 0.50 through generation " f"{above50.max() if len(above50) else 'never reached'}. " f"At generation {det_gens[-1]} the power is {pw[-1]:.3f}.") se = math.sqrt(max(pw[-1] * (1 - pw[-1]), 1e-12) / DET_REPS) P(f" Binomial standard error on that last figure: {se:.4f}.") # ----------------------------------------------------------------------- # Section 10, headline numbers # ----------------------------------------------------------------------- P("") P(banner("SECTION 10. THE NUMBERS THE ARTICLE QUOTES")) P("") r2 = [r for r in found_rows if r["k"] == 2][0] r50 = [r for r in found_rows if r["k"] == 50][0] r500 = [r for r in found_rows if r["k"] == 500][0] g2 = [g for g in grid_rows if g["k"] == 2 and g["r"] == 0.40][0] g50 = [g for g in grid_rows if g["k"] == 50 and g["r"] == 0.40][0] P(f" H retained by 2 founders : {100 * r2['ret']:.3f}% " f"+/- {100 * r2['ret_se']:.3f}% (closed form {100 * (1 - 1 / 4):.3f}%)") P(f" H retained by 50 founders : {100 * r50['ret']:.3f}% " f"+/- {100 * r50['ret_se']:.3f}% (closed form {100 * (1 - 1 / 100):.3f}%)") P(f" H retained by 500 founders : {100 * r500['ret']:.3f}% " f"+/- {100 * r500['ret_se']:.3f}% (closed form {100 * (1 - 1 / 1000):.3f}%)") P(f" alleles retained by 2 founders : {r2['alleles']:.4f} of {A_SRC:.3f} " f"per locus = {100 * r2['alleles'] / A_SRC:.2f}%") P(f" alleles retained by 50 founders : {r50['alleles']:.4f} of {A_SRC:.3f} " f"per locus = {100 * r50['alleles'] / A_SRC:.2f}%") P(f" alleles lost by 2 founders : {r2['lost']:.3f} of " f"{src_alleles_found} = {100 * r2['lost'] / src_alleles_found:.2f}%") P(f" H after growth back to K, k=2 : {100 * g2['h_end'] / H_SRC_DYN:.3f}% " f"of source, after {g2['gens']} generations at r = 0.40") P(f" H after growth back to K, k=50 : {100 * g50['h_end'] / H_SRC_DYN:.3f}% " f"of source, after {g50['gens']} generations at r = 0.40") P(f" alleles after growth back, k=2 : {g2['a_end']:.3f} per locus = " f"{100 * g2['a_end'] / A_SRC:.2f}% of source") P(f" generations for k=2 to climb half : {half[2]:,}") P(f" largest |z|, closed form : {np.abs(zs).max():.2f} over " f"{len(FOUND_SIZES)} founding sizes") P(f" largest |z|, allele loss : {np.abs(zl).max():.2f}") P(f" largest |z|, exact recursion : {np.abs(zg).max():.2f} over " f"{len(grid_rows)} cells") P(f" k = N loses nothing : max alleles lost over " f"{REPS_FINITE:,} replicates = {kN[7]:.0f}") tot_reps = (REPS_FOUND * len(FOUND_SIZES) + REPS_FINITE * len(FINITE_K) + REPS_GRID * len(grid_rows) + REPS_GRID * 3 + DET_REPS * (len(DET_FOUNDERS) + 1)) P(f" total replicates in this file : {tot_reps:,}") # ----------------------------------------------------------------------- # Section 11, figure data # ----------------------------------------------------------------------- P("") P(banner("SECTION 11. FIGURE DATA")) P("") P("Blocks below are the exact numbers plotted in the article's figures. Each") P("block begins with FIGDATA and ends with ENDDATA, one record per line.") P("") P("FIGDATA source_spectrum locus " + " ".join( f"p{i}" for i in range(K_ALLELES))) for l in range(N_LOCI_TOTAL): P(f"{l} " + " ".join(f"{v:.9f}" for v in spec[l])) P("ENDDATA") P("") P(" The block above is the complete source spectrum, every allele frequency") P(" at every locus, so that the interactive companion and any re-analysis") P(" start from the same numbers this run did.") P("") P("FIGDATA fig1_retention k retention se predicted alleles alleles_se") for r in found_rows: P(f"{r['k']} {r['ret']:.8f} {r['ret_se']:.8f} {r['pred']:.8f} " f"{r['alleles']:.6f} {r['alleles_se']:.6f}") P("ENDDATA") P("") P("FIGDATA fig2_loss k lost lost_se predicted frac_alleles_kept") for r in found_rows: P(f"{r['k']} {r['lost']:.6f} {r['lost_se']:.6f} {r['exp_lost']:.6f} " f"{r['alleles'] / A_SRC:.8f}") P("ENDDATA") P("") P("FIGDATA fig2b_perallele_k2 p pred_lost sim_lost") row = [r for r in found_rows if r["k"] == 2][0] for l in range(L_FOUND): for i in range(K_ALLELES): if spec[l, i] > 0: P(f"{spec[l, i]:.8f} {(1 - spec[l, i]) ** 4:.8f} " f"{1.0 - row['per_allele'][l, i] / REPS_FOUND:.8f}") P("ENDDATA") P("") for (k, r), tr in sorted(traj_store.items()): P(f"FIGDATA fig3_traj_k{k}_r{r:.2f} t N h h_recursion a") for i in range(len(tr["t"])): P(f"{tr['t'][i]} {tr['n'][i]} {tr['h'][i]:.8f} " f"{tr['pred'][i]:.8f} {tr['a'][i]:.6f}") P("ENDDATA") P("") P("FIGDATA fig3b_longrun k " + " ".join(str(g) for g in long_gens)) for k, vals in sorted(long_rows.items()): P(f"{k} " + " ".join(f"{v:.6f}" for v in vals)) P("ENDDATA") P("") P("FIGDATA fig4_power gen power_k5 power_k25 h_k5 h_ctrl a_k5 a_ctrl") for i, gen in enumerate(det_gens): e = det_extra[i] P(f"{gen} {det_power[5][i]:.6f} {det_power[25][i]:.6f} " f"{e[1]:.6f} {e[2]:.6f} {e[3]:.6f} {e[4]:.6f}") P("ENDDATA") P("") P("FIGDATA fig5_convergence reps k2 k10 k50") step = 25 for n in range(step, REPS_FOUND + 1, step): P(f"{n} {conv_table[2][n - 1]:.8f} {conv_table[10][n - 1]:.8f} " f"{conv_table[50][n - 1]:.8f}") P("ENDDATA") P("") P("FIGDATA grid_table k r gens h0 h0_se h_end h_se h_pred z a0 a_end") for g in grid_rows: P(f"{g['k']} {g['r']:.2f} {g['gens']} {g['h0']:.8f} {g['h0_se']:.8f} " f"{g['h_end']:.8f} {g['h_se']:.8f} {g['h_pred']:.8f} {g['z']:.4f} " f"{g['a0']:.6f} {g['a_end']:.6f}") P("ENDDATA") P("") P("FIGDATA spectrum_shift class_lo class_hi source founders_k5 recovered_k5") # spectrum of founders and of the recovered population, k = 5, r = 0.40 rng_sp = np.random.default_rng(streams[6]) sizes_sp = size_path_geometric(5, 0.40, K_CAP, GRID_TMIN, GRID_TMAX, GRID_SETTLE) reps_sp = 2000 p_sp_tile = np.repeat(p_det[None, :, :], reps_sp, axis=0).reshape(-1, K_ALLELES) cnt_sp = rng_sp.multinomial(2 * 5, p_sp_tile).astype(np.int32) f_found = (cnt_sp / cnt_sp.sum(axis=1, keepdims=True)).ravel() cnt_sp2 = evolve(rng_sp, cnt_sp, sizes_sp, MU, K_ALLELES) f_rec = (cnt_sp2 / cnt_sp2.sum(axis=1, keepdims=True)).ravel() src_flat2 = p_det.ravel() hs_src = np.histogram(src_flat2[src_flat2 > 0], bins=edges)[0] / float(L_DET) hs_f = np.histogram(f_found[f_found > 0], bins=edges)[0] / float(reps_sp * L_DET) hs_r = np.histogram(f_rec[f_rec > 0], bins=edges)[0] / float(reps_sp * L_DET) for i in range(len(edges) - 1): P(f"{edges[i]:.4f} {min(edges[i + 1], 1.0):.4f} {hs_src[i]:.6f} " f"{hs_f[i]:.6f} {hs_r[i]:.6f}") P("ENDDATA") P("") P(f" Spectrum note: the source carries {hs_src[0]:.3f} alleles per locus below 1%") P(f" frequency. Five founders keep {hs_f[0]:.3f} of them, and after " f"{len(sizes_sp)} generations") P(f" back at carrying capacity the population has {hs_r[0]:.3f} in that class, " "nearly all of") P(" them new mutations rather than survivors.") P("") P(rule()) P(f"total wall clock: {time.time() - t_start:.1f} s") P(rule()) sys.stdout.write("\n".join(out) + "\n") if __name__ == "__main__": main()