"""Part 1: are 'singular' evolutionary innovations older than repeated ones, and what does singularity tell us about hard steps? Inputs : data/vermeij2006_recoded_2026.csv, data/innovations_2026.csv (see build_innovations.py) Outputs: data/part1_tests.csv, data/part1_recurrence_rates.csv, data/part1_summary.json, figures/fig1_singular_vs_repeated_age.png, figures/fig2_recurrence_rates.png """ import csv import json import os import numpy as np from scipy import stats import plotstyle HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.join(HERE, "..", "data") FIG = os.path.join(HERE, "..", "figures") RNG = np.random.default_rng(20261011) NPERM = 20000 # ------------------------------------------------------------------ statistics helpers def firth_logistic(x, y, maxit=200, tol=1e-10): """Firth (1993) bias-reduced logistic regression of y on [1, x]. Returns beta (2,), covariance (2,2), penalised log-likelihood.""" X = np.column_stack([np.ones_like(x), x]) b = np.zeros(2) for _ in range(maxit): eta = X @ b p = 1.0 / (1.0 + np.exp(-eta)) W = p * (1 - p) XtWX = X.T @ (X * W[:, None]) inv = np.linalg.inv(XtWX) h = np.einsum("ij,jk,ik->i", X * np.sqrt(W)[:, None], inv, X * np.sqrt(W)[:, None]) U = X.T @ (y - p + h * (0.5 - p)) step = inv @ U step = np.clip(step, -5, 5) b = b + step if np.max(np.abs(step)) < tol: break eta = X @ b p = 1.0 / (1.0 + np.exp(-eta)) W = p * (1 - p) XtWX = X.T @ (X * W[:, None]) ll = np.sum(y * np.log(p + 1e-300) + (1 - y) * np.log(1 - p + 1e-300)) + 0.5 * np.linalg.slogdet(XtWX)[1] return b, np.linalg.inv(XtWX), ll def firth_profile_ci(x, y, level=3.841): """Profile penalised-likelihood 95 % interval for the slope of the Firth logistic regression (the interval recommended for Firth fits; Heinze & Schemper 2002). For each fixed slope the intercept is maximised numerically.""" from scipy import optimize X1 = np.ones_like(x) def pll(b0, b1): eta = b0 + b1 * x p = 1.0 / (1.0 + np.exp(-eta)) W = p * (1 - p) a, b, c = np.sum(W * X1), np.sum(W * x), np.sum(W * x * x) det = a * c - b * b if det <= 0: return -1e300 return float(np.sum(y * np.log(p + 1e-300) + (1 - y) * np.log(1 - p + 1e-300)) + 0.5 * np.log(det)) def prof(b1): r = optimize.minimize_scalar(lambda b0: -pll(b0, b1), bounds=(-40, 40), method="bounded", options=dict(xatol=1e-8)) return -r.fun bhat, _, _ = firth_logistic(x, y) top = prof(bhat[1]) def g(b1): return 2 * (top - prof(b1)) - level lo = hi = np.nan try: lo = optimize.brentq(g, bhat[1] - 30, bhat[1]) except ValueError: pass try: hi = optimize.brentq(g, bhat[1], bhat[1] + 30) except ValueError: pass return lo, hi def tests(age, sing, label): """All tests for one coding. age in Ma (array), sing boolean array (True = singular).""" age = np.asarray(age, float) sing = np.asarray(sing, bool) la = np.log10(age) out = dict(coding=label, n=int(len(age)), n_singular=int(sing.sum()), n_repeated=int((~sing).sum())) if sing.sum() < 2 or (~sing).sum() < 2: return out out["median_age_singular_Ma"] = float(np.median(age[sing])) out["median_age_repeated_Ma"] = float(np.median(age[~sing])) for cut, nm in ((1000.0, "1Ga"), (543.0, "543Ma")): a = int((sing & (age >= cut)).sum()); b = int((sing & (age < cut)).sum()) c = int((~sing & (age >= cut)).sum()); d = int((~sing & (age < cut)).sum()) odds, p = stats.fisher_exact([[a, b], [c, d]], alternative="two-sided") out[f"fisher_{nm}_table"] = f"{a}/{a+b} vs {c}/{c+d}" out[f"fisher_{nm}_p"] = float(p) u = stats.mannwhitneyu(la[sing], la[~sing], alternative="two-sided") out["mannwhitney_p"] = float(u.pvalue) # permutation test on the difference of mean log10 age (two-sided) obs = la[sing].mean() - la[~sing].mean() k = int(sing.sum()) cnt = 0 for _ in range(NPERM): perm = RNG.permutation(la) if abs(perm[:k].mean() - perm[k:].mean()) >= abs(obs) - 1e-12: cnt += 1 out["diff_mean_log10age"] = float(obs) out["perm_p"] = (cnt + 1) / (NPERM + 1) b, cov, ll = firth_logistic(la, sing.astype(float)) out["firth_slope_per_dex"] = float(b[1]) out["firth_slope_se"] = float(np.sqrt(cov[1, 1])) out["odds_ratio_per_dex"] = float(np.exp(b[1])) out["odds_ratio_lo95"] = float(np.exp(b[1] - 1.96 * np.sqrt(cov[1, 1]))) out["odds_ratio_hi95"] = float(np.exp(b[1] + 1.96 * np.sqrt(cov[1, 1]))) plo, phi = firth_profile_ci(la, sing.astype(float)) out["odds_ratio_profile_lo95"] = float(np.exp(plo)) if np.isfinite(plo) else None out["odds_ratio_profile_hi95"] = float(np.exp(phi)) if np.isfinite(phi) else None # permutation p-value for the slope cnt = 0 nper = 2000 for _ in range(nper): bp, _, _ = firth_logistic(la, RNG.permutation(sing.astype(float))) if abs(bp[1]) >= abs(b[1]) - 1e-12: cnt += 1 out["firth_slope_perm_p"] = (cnt + 1) / (nper + 1) # leave-one-out range of the slope loo = [] for i in range(len(la)): m = np.ones(len(la), bool); m[i] = False if sing[m].sum() < 2 or (~sing[m]).sum() < 2: continue loo.append(firth_logistic(la[m], sing[m].astype(float))[0][1]) out["loo_slope_min"] = float(np.min(loo)); out["loo_slope_max"] = float(np.max(loo)) out["_beta"] = [float(b[0]), float(b[1])] return out def load(): v = list(csv.DictReader(open(os.path.join(DATA, "vermeij2006_recoded_2026.csv"), encoding="utf8"))) e = list(csv.DictReader(open(os.path.join(DATA, "innovations_2026.csv"), encoding="utf8"))) return v, e def main(): plt = plotstyle.apply() v, e = load() res = [] # ---------- A. Vermeij as published, and our 2026 re-codings ageP = np.array([float(r["age_published_Ma"]) for r in v]) age26 = np.array([float(r["age_2026_Ma"]) for r in v]) sP = np.array([r["status_published"] == "singular" for r in v]) sS = np.array([r["status_2026_strict"] == "singular" for r in v]) sF = np.array([r["status_2026_functional"] == "singular" for r in v]) res.append(tests(ageP, sP, "Vermeij 2006 as published (n=77)")) res.append(tests(age26, sS, "Vermeij list, 2026 strict re-coding")) res.append(tests(age26, sF, "Vermeij list, 2026 functional re-coding")) # subsets of the published list ph = ageP < 543 res.append(tests(ageP[ph], sP[ph], "Vermeij as published, Phanerozoic only (<543 Ma)")) ph2 = (ageP < 543) & (ageP > 5) res.append(tests(ageP[ph2], sP[ph2], "Vermeij as published, Phanerozoic, excluding <5 Ma (language)")) nl = ageP < 3000 res.append(tests(ageP[nl], sP[nl], "Vermeij as published, excluding LUCA-era items (>=3.5 Ga)")) pc = ageP >= 543 res.append(tests(ageP[pc], sP[pc], "Vermeij as published, Precambrian only (>=543 Ma)")) # co-inherited stem-lineage traits counted once (referee M2/M3). A cluster is replaced by one row at # the cluster's median age; rows outside any cluster are unchanged. def collapse(col): keep_age, keep_s, seen = [], [], {} for r, a, st in zip(v, ageP, sP): c = r[col] if not c: keep_age.append(a); keep_s.append(st) else: seen.setdefault(c, []).append(a) for c, ages in seen.items(): keep_age.append(float(np.median(ages))); keep_s.append(True) return np.array(keep_age), np.array(keep_s) def only_bodyplan(): keep_age, keep_s, bp = [], [], [] for r, a, st in zip(v, ageP, sP): if r["cluster_major"] == "animal body plan": bp.append(a) else: keep_age.append(a); keep_s.append(st) keep_age.append(float(np.median(bp))); keep_s.append(True) return np.array(keep_age), np.array(keep_s) collapsed = {} for lab, (ca, cs) in (("animal body-plan rows (7) counted as one event", only_bodyplan()), ("stem-lineage clusters counted once (LUCA, LECA, Eumetazoa stem, Bilateria stem)", collapse("cluster_stem")), ("LUCA (4), LECA (2) and animal body plan (7) each counted as one event", collapse("cluster_major"))): t_ = tests(ca, cs, "Vermeij as published, " + lab) res.append(t_) collapsed[lab] = dict(n=t_["n"], n_singular=t_["n_singular"], fisher_543Ma_table=t_["fisher_543Ma_table"], fisher_543Ma_p=t_["fisher_543Ma_p"], fisher_1Ga_table=t_["fisher_1Ga_table"], fisher_1Ga_p=t_["fisher_1Ga_p"], odds_ratio_per_dex=t_["odds_ratio_per_dex"], odds_ratio_profile_ci95=[t_["odds_ratio_profile_lo95"], t_["odds_ratio_profile_hi95"]], firth_slope_perm_p=t_["firth_slope_perm_p"]) ph26 = age26 < 543 res.append(tests(age26[ph26], sS[ph26], "Vermeij 2026 strict, Phanerozoic only")) res.append(tests(age26[ph26], sF[ph26], "Vermeij 2026 functional, Phanerozoic only")) # ---------- B. extended list of general innovations ee = [r for r in e if r["age_best_Ma"] not in ("", "None")] ab = np.array([float(r["age_best_Ma"]) for r in ee]) amin = np.array([float(r["age_min_Ma"]) for r in ee]) amax = np.array([float(r["age_max_Ma"]) for r in ee]) nb = np.array([int(r["n_best"]) for r in ee]) nmin = np.array([int(r["n_min"]) for r in ee]) nmax = np.array([int(r["n_max"]) for r in ee]) eS = nb == 1 eF = nmax == 1 res.append(tests(ab, eS, f"Extended list strict (n={len(ee)}), best ages")) res.append(tests(ab, eF, "Extended list functional, best ages")) res.append(tests(amin, eS, "Extended list strict, youngest ages")) res.append(tests(amax, eS, "Extended list strict, oldest ages")) old = ab > 5 res.append(tests(ab[old], eS[old], "Extended list strict, excluding innovations <5 Ma old")) res.append(tests(ab[old], eF[old], "Extended list functional, excluding innovations <5 Ma old")) ph = (ab < 543) & (ab > 5) res.append(tests(ab[ph], eS[ph], "Extended list strict, Phanerozoic (5-543 Ma)")) # the four singular rows younger than 3 Ma are one lineage (Homo): count them once (median age) homo = np.array([r["lineage_cluster"] == "Homo" for r in ee]) ab1 = np.concatenate([ab[~homo], [np.median(ab[homo])]]) eS1 = np.concatenate([eS[~homo], [True]]) res.append(tests(ab1, eS1, "Extended list strict, the four Homo rows counted as one event")) # Monte-Carlo over coding and dating uncertainty in the extended list slopes, pvals = [], [] for _ in range(2000): a = 10 ** RNG.uniform(np.log10(amin), np.log10(amax)) # contested items (n_best==1 but n_max>1, or n_min==1 but n_best>1) are coded at random s = np.where(nmax == 1, True, np.where(nmin > 1, False, RNG.random(len(nb)) < 0.5)) if s.sum() < 2 or (~s).sum() < 2: continue b, cov, _ = firth_logistic(np.log10(a), s.astype(float)) slopes.append(b[1]) pvals.append(stats.mannwhitneyu(np.log10(a[s]), np.log10(a[~s])).pvalue) mc = dict(n_draws=len(slopes), slope_median=float(np.median(slopes)), slope_p05=float(np.percentile(slopes, 5)), slope_p95=float(np.percentile(slopes, 95)), frac_slope_positive=float(np.mean(np.array(slopes) > 0)), frac_mannwhitney_p_below_0p05=float(np.mean(np.array(pvals) < 0.05))) # same Monte-Carlo excluding innovations younger than 5 Ma slopes2 = [] keep = ab > 5 for _ in range(2000): a = 10 ** RNG.uniform(np.log10(amin[keep]), np.log10(amax[keep])) s = np.where(nmax[keep] == 1, True, np.where(nmin[keep] > 1, False, RNG.random(keep.sum()) < 0.5)) if s.sum() < 2 or (~s).sum() < 2: continue slopes2.append(firth_logistic(np.log10(a), s.astype(float))[0][1]) mc["excl_young_slope_median"] = float(np.median(slopes2)) mc["excl_young_slope_p05"] = float(np.percentile(slopes2, 5)) mc["excl_young_slope_p95"] = float(np.percentile(slopes2, 95)) mc["excl_young_frac_slope_positive"] = float(np.mean(np.array(slopes2) > 0)) # ---------- C. information-loss diagnostics # (i) among repeated innovations, does the number of documented origins decline with age? rep = nb >= 2 rho, p_rho = stats.spearmanr(np.log10(ab[rep]), np.log10(nb[rep])) cnt = 0 for _ in range(NPERM): r_, _ = stats.spearmanr(RNG.permutation(np.log10(ab[rep])), np.log10(nb[rep])) if abs(r_) >= abs(rho) - 1e-12: cnt += 1 info = dict(spearman_rho_norigins_vs_age=float(rho), spearman_p=float(p_rho), spearman_perm_p=(cnt + 1) / (NPERM + 1), n_repeated=int(rep.sum())) # (ii) gradient vs step: compare a log-age gradient with a step at 1 Ga (published list) la = np.log10(ageP) y = sP.astype(float) def ll_step(cut): ll = 0.0 for m in (ageP >= cut, ageP < cut): k, n = y[m].sum(), m.sum() p = (k + 0.5) / (n + 1.0) ll += k * np.log(p) + (n - k) * np.log(1 - p) return ll b, cov, _ = firth_logistic(la, y) p = 1 / (1 + np.exp(-(b[0] + b[1] * la))) ll_grad = float(np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))) info["loglik_gradient_logage"] = ll_grad info["loglik_step_1Ga"] = float(ll_step(1000.0)) info["loglik_step_543Ma"] = float(ll_step(543.0)) # scan of the cut position (referee M3): for which cuts does a step fit better than the gradient? ua = np.unique(ageP) cuts = np.sqrt(ua[:-1] * ua[1:]) # geometric midpoints between adjacent distinct ages scan = [dict(cut_Ma=float(c), between_Ma=[float(a_), float(b_)], loglik_step=float(ll_step(c)), singular_older=f"{int(y[ageP >= c].sum())}/{int((ageP >= c).sum())}", singular_younger=f"{int(y[ageP < c].sum())}/{int((ageP < c).sum())}") for c, a_, b_ in zip(cuts, ua[:-1], ua[1:])] info["step_cut_scan"] = scan def runs(thr): """Contiguous runs of cut positions where the step log-likelihood exceeds the gradient's by more than thr.""" out, cur = [], None for d in scan: if d["loglik_step"] - ll_grad > thr: cur = [d["between_Ma"][0], d["between_Ma"][1]] if cur is None else [cur[0], d["between_Ma"][1]] elif cur is not None: out.append(cur); cur = None if cur is not None: out.append(cur) return out info["step_minus_gradient_loglik_exceeds"] = {"0": runs(0.0), "1": runs(1.0), "2": runs(2.0), "3": runs(3.0)} info["step_best_cut"] = max(scan, key=lambda d: d["loglik_step"]) info["step_note"] = ("the cut is a third fitted parameter, so a fair comparison charges the step about one log unit (AIC) " "or more; the intervals above are the cut positions (in Ma, between adjacent dated rows) for which " "the step log-likelihood exceeds the two-parameter gradient's by more than 0, 1, 2 and 3 units") def step_vs_gradient(age_, y_): age_ = np.asarray(age_, float); y_ = np.asarray(y_, float) b_, _, _ = firth_logistic(np.log10(age_), y_) p_ = 1 / (1 + np.exp(-(b_[0] + b_[1] * np.log10(age_)))) lg = float(np.sum(y_ * np.log(p_) + (1 - y_) * np.log(1 - p_))) ua_ = np.unique(age_) best = None for c in np.sqrt(ua_[:-1] * ua_[1:]): l = 0.0 for m in (age_ >= c, age_ < c): k, n = y_[m].sum(), m.sum() q = (k + 0.5) / (n + 1.0) l += k * np.log(q) + (n - k) * np.log(1 - q) if best is None or l > best[1]: best = (float(c), float(l)) return dict(loglik_gradient=lg, loglik_best_step=best[1], best_cut_Ma=best[0], difference=best[1] - lg) info["step_vs_gradient_after_collapsing"] = { "animal body-plan rows (7) counted as one event": step_vs_gradient(*only_bodyplan()), "stem-lineage clusters counted once": step_vs_gradient(*collapse("cluster_stem")), "LUCA, LECA and animal body plan each counted once": step_vs_gradient(*collapse("cluster_major"))} info["collapsed_coinherited"] = collapsed info["vermeij_text_vs_table"] = dict( text="14 of 23 singular (61%) and 9 of 55 repeated (16%) before 543 Ma; 8 (35%) and 4 (7.3%) before 1 Ga; P < 0.02 and P < 0.01", table_as_printed="15 of 23 singular and 9 of 54 repeated dated 550 Ma or older; 8 and 4 at 1 Ga or older", note="the printed tables give one more singular row before 543 Ma than the text, and one fewer repeated row in total") # proportion singular by era (published list) eras = [("<66 Ma", 0, 66), ("66-252 Ma", 66, 252), ("252-543 Ma", 252, 543), ("543-1000 Ma", 543, 1000), (">=1000 Ma", 1000, 1e5)] era_tab = [] for nm, lo, hi in eras: m = (ageP >= lo) & (ageP < hi) m26 = (age26 >= lo) & (age26 < hi) era_tab.append(dict(era=nm, n_published=int(m.sum()), singular_published=int(sP[m].sum()), n_2026=int(m26.sum()), singular_2026_strict=int(sS[m26].sum()), singular_2026_functional=int(sF[m26].sum()))) info["era_table"] = era_tab # (iii) ILLUSTRATION ONLY (demoted in revision, referee M2): the ratio of the Phanerozoic singular # fraction to the Precambrian one. It could be read as P(genuinely unique | appears singular, # Precambrian) only if (a) Phanerozoic and Precambrian rows were exchangeable kinds of innovation, # (b) 'singular' in the Phanerozoic meant 'improbable', and (c) rows were independent. The list # violates all three (minor clade-defined Phanerozoic traits against major Precambrian transitions; # co-inherited stem-lineage rows), so the number is NOT an estimate; the true value could lie # anywhere from 0 to 1. The binomial interval below is sampling error only. def base_rate(k1, n1, k2, n2): q = RNG.beta(k1 + 1, n1 - k1 + 1, 200000) f = RNG.beta(k2 + 1, n2 - k2 + 1, 200000) pi = np.minimum(q / f, 1.0) return dict(phanerozoic_singular=int(k1), phanerozoic_n=int(n1), precambrian_singular=int(k2), precambrian_n=int(n2), q_phan=float(k1 / n1), f_pre=float(k2 / n2), ratio=float(min((k1 / n1) / (k2 / n2), 1.0)) if k2 else None, ratio_posterior_median=float(np.median(pi)), binomial_only_interval95=[float(np.percentile(pi, 2.5)), float(np.percentile(pi, 97.5))]) br = {} m1 = (ageP > 5) & (ageP < 543) m2 = ageP >= 543 br["Vermeij as published"] = base_rate(sP[m1].sum(), m1.sum(), sP[m2].sum(), m2.sum()) m1 = (age26 > 5) & (age26 < 543) m2 = age26 >= 543 br["Vermeij, 2026 strict re-coding"] = base_rate(sS[m1].sum(), m1.sum(), sS[m2].sum(), m2.sum()) br["Vermeij, 2026 functional re-coding"] = base_rate(sF[m1].sum(), m1.sum(), sF[m2].sum(), m2.sum()) ed = [r for r in e if r["age_best_Ma"] not in ("", "None")] ea = np.array([float(r["age_best_Ma"]) for r in ed]) es_ = np.array([int(r["n_best"]) == 1 for r in ed]) m1 = (ea > 5) & (ea < 543) m2 = ea >= 543 br["Extended 2026 list (strict)"] = base_rate(es_[m1].sum(), m1.sum(), es_[m2].sum(), m2.sum()) # range of the ratio across cut positions and collapses of co-inherited rows (published list) variants = {} for cut in (543.0, 610.0, 1000.0): for lab, (ca, cs) in (("all rows", (ageP, sP)), ("body-plan rows as one", only_bodyplan()), ("stem clusters once", collapse("cluster_stem")), ("major clusters once", collapse("cluster_major"))): for excl_young in (True, False): m1 = (ca < cut) & ((ca > 5) if excl_young else True) m2 = ca >= cut if cs[m2].sum() == 0: continue variants[f"cut {cut:.0f} Ma; {lab}; {'excluding' if excl_young else 'including'} language (1 Ma)"] = float( min((cs[m1].sum() / m1.sum()) / (cs[m2].sum() / m2.sum()), 1.0)) info["singular_fraction_ratio_illustration"] = dict( by_coding=br, variants=variants, variant_range=[min(variants.values()), max(variants.values())], warning="illustration under an untested exchangeability assumption; not an estimate of the share of " "ancient singular innovations that are genuinely improbable, which could be anywhere from 0 to 1") # ---------- D. recurrence rates (opportunity-corrected) rr = [] for r in e: if r["age_best_Ma"] in ("", "None"): continue a, a_lo, a_hi = float(r["age_best_Ma"]), float(r["age_min_Ma"]), float(r["age_max_Ma"]) n, n_lo, n_hi = int(r["n_best"]), int(r["n_min"]), int(r["n_max"]) d = dict(id=r["id"], innovation=r["innovation"], set=r["set"], cls=r["class"], n_best=n, age_best_Ma=a) if n == 1: # no recurrence in a Myr: 95% upper limit on a Poisson recurrence rate = 3.0/a d["rate_per_Gyr"] = None d["rate_upper95_per_Gyr"] = 2.996 / a * 1000.0 d["rate_upper95_per_Gyr_if_oldest"] = 2.996 / a_hi * 1000.0 d["rate_upper95_per_Gyr_if_youngest"] = 2.996 / a_lo * 1000.0 # likelihood of zero recurrences if the true recurrence rate were 1 per Gyr / 0.2 per Gyr d["P0_rate1perGyr"] = float(np.exp(-a / 1000.0)) d["P0_rate0p2perGyr"] = float(np.exp(-0.2 * a / 1000.0)) # the limit depends on age only: it is the same whatever the true rate, and lies above the # hard-step regime (< ~0.2 per Gyr) for anything younger than 15 Gyr, so it cannot test it d["limit_above_hardstep_regime"] = int(2.996 / a * 1000.0 > 0.2) else: d["rate_per_Gyr"] = (n - 1) / a * 1000.0 d["rate_lo_per_Gyr"] = max(n_lo - 1, 0.5) / a_hi * 1000.0 d["rate_hi_per_Gyr"] = (n_hi - 1) / a_lo * 1000.0 rr.append(d) # ---------- write tables keys = [] for r in res: for k in r: if k not in keys and not k.startswith("_"): keys.append(k) with open(os.path.join(DATA, "part1_tests.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=keys, extrasaction="ignore") w.writeheader() for r in res: w.writerow(r) rkeys = ["id", "innovation", "set", "cls", "n_best", "age_best_Ma", "rate_per_Gyr", "rate_lo_per_Gyr", "rate_hi_per_Gyr", "rate_upper95_per_Gyr", "rate_upper95_per_Gyr_if_oldest", "rate_upper95_per_Gyr_if_youngest", "P0_rate1perGyr", "P0_rate0p2perGyr", "limit_above_hardstep_regime"] with open(os.path.join(DATA, "part1_recurrence_rates.csv"), "w", newline="") as f: w = csv.DictWriter(f, fieldnames=rkeys, extrasaction="ignore") w.writeheader() w.writerows(rr) hs = [r for r in e if r["set"] == "hardstep"] summary = dict(tests=[{k: v_ for k, v_ in r.items() if not k.startswith("_")} for r in res], monte_carlo_extended=mc, information_loss=info, hardstep_class_counts={c: sum(1 for r in hs if r["class"] == c) for c in "ABC"}, hardstep_classes={r["id"]: r["class"] for r in hs}) json.dump(summary, open(os.path.join(DATA, "part1_summary.json"), "w"), indent=1) # ---------- Figure 1: ages of singular vs repeated innovations fig, axes = plt.subplots(1, 2, figsize=(11, 4.4), gridspec_kw={"width_ratios": [1.25, 1]}) ax = axes[0] jit = RNG.uniform(-0.16, 0.16, len(ageP)) clustered = np.array([bool(r["cluster_major"]) for r in v]) for flag, yv, col, mk, lab in ((True, 1, plotstyle.ORANGE, "o", "claimed unique (23)"), (False, 0, plotstyle.BLUE, "s", "repeated (54)")): m = (sP == flag) & ~clustered ax.scatter(ageP[m], yv + jit[m], s=34, color=col, marker=mk, edgecolor=plotstyle.SURFACE, linewidth=0.8, label=lab, zorder=3) m = clustered ax.scatter(ageP[m], 1 + jit[m], s=34, facecolors=plotstyle.SURFACE, edgecolors=plotstyle.ORANGE, marker="o", linewidth=1.4, zorder=3) ax.text(300, 0.66, "hollow circles: 13 co-inherited rows\n(4 LUCA, 2 LECA, 7 animal body plan)", fontsize=8, color=plotstyle.INK2, ha="left", va="center") r0 = res[0] xs = np.logspace(0, np.log10(4500), 200) b0, b1 = r0["_beta"] ax.plot(xs, 1 / (1 + np.exp(-(b0 + b1 * np.log10(xs)))), color=plotstyle.INK2, lw=1.6, zorder=2) ax.text(40, 0.30, "fitted P(unique | age)", color=plotstyle.INK2, fontsize=8.5, ha="left") ax.axvline(543, color=plotstyle.GREY, lw=1, ls="--", zorder=1) ax.text(620, 1.34, "Precambrian", fontsize=8, color=plotstyle.INK2, ha="right", va="top") ax.text(470, 1.34, "Phanerozoic", fontsize=8, color=plotstyle.INK2, ha="left", va="top") ax.set_xscale("log") ax.set_xlim(4500, 0.8) ax.set_ylim(-0.35, 1.38) ax.set_yticks([0, 1]) ax.set_yticklabels(["repeated\n(54)", "claimed\nunique (23)"]) ax.set_xlabel("Age of first origin (Ma, log scale; older to the left)") ax.set_title("Vermeij (2006) innovations, as published:\nunique ones are older, mostly through co-inherited early rows") ax.grid(axis="y", visible=False) ax = axes[1] labs = [d["era"] for d in era_tab] xpos = np.arange(len(labs)) fp = [d["singular_published"] / d["n_published"] if d["n_published"] else np.nan for d in era_tab] f26 = [d["singular_2026_strict"] / d["n_2026"] if d["n_2026"] else np.nan for d in era_tab] ff = [d["singular_2026_functional"] / d["n_2026"] if d["n_2026"] else np.nan for d in era_tab] w = 0.26 ax.bar(xpos - w, fp, w * 0.92, color=plotstyle.ORANGE, label="as published") ax.bar(xpos, f26, w * 0.92, color=plotstyle.BLUE, label="2026 strict re-coding") ax.bar(xpos + w, ff, w * 0.92, color=plotstyle.AQUA, label="2026 functional re-coding") for i, d in enumerate(era_tab): ax.text(i - w, (fp[i] if fp[i] == fp[i] else 0) + 0.02, f"{d['singular_published']}/{d['n_published']}", ha="center", fontsize=7.5, color=plotstyle.INK2) ax.text(i + w, (ff[i] if ff[i] == ff[i] else 0) + 0.02, f"{d['singular_2026_functional']}/{d['n_2026']}", ha="center", fontsize=7.5, color=plotstyle.INK2) ax.set_xticks(xpos) ax.set_xticklabels(labs, fontsize=8) ax.set_ylim(0, 1.12) ax.set_ylabel("Fraction of innovations that are unique") ax.set_title("Fraction unique, by era of first origin") ax.legend(loc="upper left") ax.grid(axis="x", visible=False) fig.tight_layout() fig.savefig(os.path.join(FIG, "fig1_singular_vs_repeated_age.png")) plt.close(fig) # ---------- Figure 2: recurrence rates fig, ax = plt.subplots(figsize=(10.2, 6.6)) # Selective direct labels: every proposed hard step and a few comparison innovations are labelled; # the remaining comparison points are listed in data/part1_recurrence_rates.csv. LAB = {"language": (6, 3, "left"), "hominin_brain": (6, 3, "left"), "fire": (6, 3, "left"), "cum_culture": (6, 3, "left"), "c4": (6, 3, "left"), "viviparity": (6, 3, "left"), "agriculture": (6, 3, "left"), "big_brain": (-6, 11, "right"), "nitroplast": (6, 3, "left"), "land": (-7, -2, "right"), "abiogenesis": (9, -3, "left"), "phototrophy": (-7, -6, "right"), "dna_replication": (7, -11, "left"), "simple_multicell": (8, 3, "left")} # labels placed in free space with a leader line: id -> (x_text, y_text, ha) TXT = {"complex_multicell": (1000, 40, "left"), "animal_multicell": (200, 3.6, "right"), "neurons": (200, 2.4, "right"), "eukaryogenesis": (200, 1.6, "right"), "oxyphot": (200, 1.07, "right"), "plastid": (200, 0.7, "right"), "n_fixation": (8000, 2.6, "left")} NAME = {"eukaryogenesis": "eukaryogenesis = mitochondrion = meiosis", "mitochondrion": "", "meiosis": "", "oxyphot": "oxygenic photosynthesis", "cum_culture": "cumulative culture", "c4": "C4 photosynthesis", "n_fixation": "N fixation", "dna_replication": "DNA replication", "animal_multicell": "animal multicellularity", "complex_multicell": "complex multicellularity", "simple_multicell": "simple multicellularity", "big_brain": "large brains", "land": "life on land", "hominin_brain": "hominin-level brain"} for d in rr: a = d["age_best_Ma"] hard = d["set"] == "hardstep" col = plotstyle.ORANGE if hard else plotstyle.BLUE if d["n_best"] == 1: yv = d["rate_upper95_per_Gyr"] ax.scatter([a], [yv], marker="v", s=64, color=col, edgecolor=plotstyle.SURFACE, linewidth=0.8, zorder=4) ax.annotate("", xy=(a, yv / 3.2), xytext=(a, yv), arrowprops=dict(arrowstyle="-", color=col, lw=1.2)) else: yv = d["rate_per_Gyr"] ax.errorbar([a], [yv], yerr=[[yv - d["rate_lo_per_Gyr"]], [max(d["rate_hi_per_Gyr"] - yv, 0)]], fmt="o", ms=6, color=col, mec=plotstyle.SURFACE, mew=0.8, elinewidth=1, zorder=3) nm = NAME.get(d["id"], d["id"].replace("_", " ")) if d["id"] in TXT: tx, ty, ha = TXT[d["id"]] ax.annotate(nm, (a, yv), xytext=(tx, ty), textcoords="data", fontsize=7.6, color=plotstyle.INK2, ha=ha, va="center", arrowprops=dict(arrowstyle="-", color=plotstyle.GREY, lw=0.6, shrinkB=3)) elif d["id"] in LAB: dx, dy, ha = LAB[d["id"]] ax.annotate(nm, (a, yv), xytext=(dx, dy), textcoords="offset points", fontsize=7.6, color=plotstyle.INK2, ha=ha) xs = np.logspace(-1, 3.7, 50) ax.plot(xs, 2.996 / xs * 1000, color=plotstyle.GREY, lw=1, ls="--", zorder=1) ax.annotate("dashed line: 95% upper limit for a\nunique innovation = 3 / age", xy=(15, 2.996 / 15 * 1000), xytext=(9, 3.0e4), fontsize=8, color=plotstyle.INK2, arrowprops=dict(arrowstyle="-", color=plotstyle.INK2, lw=0.7)) ax.axhline(0.2, color=plotstyle.RED, lw=1, ls=":", zorder=1) ax.text(0.22, 0.24, "hard-step regime: rate below ~1 per 5 Gyr", fontsize=8, color=plotstyle.INK2) ax.text(0.22, 0.075, "every upper limit lies above this line by construction, so the limits cannot test hard steps", fontsize=8, color=plotstyle.INK2) ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlim(0.2, 30000); ax.set_ylim(0.05, 3e5) ax.set_xticks([1, 10, 100, 1000, 10000]); ax.set_xticklabels(["1", "10", "100", "1,000", "10,000"]) ax.set_xlabel("Age of first origin (Ma)") ax.set_ylabel("Recurrence rate after first origin (independent origins per Gyr)") ax.set_title("How often did each innovation recur after its first origin?\nOld unique innovations cannot be told apart from old repeated ones") from matplotlib.lines import Line2D ax.legend(handles=[Line2D([], [], marker="s", ls="", color=plotstyle.ORANGE, label="proposed hard step"), Line2D([], [], marker="s", ls="", color=plotstyle.BLUE, label="comparison innovation"), Line2D([], [], marker="o", ls="", color=plotstyle.INK2, label="repeated: (origins - 1) / age"), Line2D([], [], marker="v", ls="", color=plotstyle.INK2, label="unique: 95% upper limit")], loc="upper right") fig.tight_layout() fig.savefig(os.path.join(FIG, "fig2_recurrence_rates.png")) plt.close(fig) return summary if __name__ == "__main__": s = main() for t in s["tests"]: print(t["coding"]) print(" ", {k: (round(v, 4) if isinstance(v, float) else v) for k, v in t.items() if k != "coding"}) print(json.dumps(s["monte_carlo_extended"], indent=1)) print(json.dumps(s["information_loss"], indent=1)) print(s["hardstep_class_counts"], s["hardstep_classes"])