"""Transfer step, rewritten after the referee report. What this module does NOT do any more: recommend a prior for the project model's expansion probability. The referee showed that (a) the lower tail of the first-pass 'analogue-informed prior' was set by an arbitrary floor on a judgement- based feasibility discount, (b) the reference classes behind the three behavioural factors were chosen after the fact, and (c) the quantity the analogues measure is not the model's parameter. What it does: 1. counts the maritime census under a stated rule, at episode and at lineage level, with a sensitivity test for capable cultures that are missing from the list; 2. forms an ILLUSTRATIVE product (launch x keeps going x persists), conditional on feasibility and valid only for the first few relaunch generations in a closed domain, under four explicit scenarios; 3. shows how the first-pass 'prior' moves with the feasibility floor (1e-2, 1e-4, 1e-6); 4. re-runs the project model with several ASSUMED expansion priors, as a sensitivity test only; 5. quantifies the mismatch between the analogue quantity and the model's parameter (targets, generations, duration) and how many observed wave-generations it would take to bound the per-generation stopping hazard usefully. JUDGEMENT marks every number that is a choice, not data. """ import os import sys import numpy as np import pandas as pd from common import CUR, DATA, HERE, SEED, dump_json C_LIGHT = 2.998e8 LY_KM = 9.4607e12 def cultures(): t = pd.read_csv(os.path.join(CUR, "maritime_cultures.csv")) c = t[t.counted == 1] out = dict(rule="culture-period with documented open-sea voyages and uninhabited, later-habitable land within about twice its " "longest documented crossing; rule written after the list existed; list not exhaustive", rows_listed=int(len(t)), rows_counted=int(len(c)), rows_not_counted={r.culture: r.note.split(":")[0] for _, r in t[t.counted == 0].iterrows()}, grades_counted=c.grade.value_counts().to_dict()) for col in ("outcome_strict", "outcome_lenient"): v = c[col].dropna() k, n = int(v.sum()), int(len(v)) out[col] = dict(k=k, n=n, fraction=k / n) # lineage level: a lineage counts as an expander if any of its episodes launched a multi-landmass wave lin = c.groupby("lineage").outcome_strict.max() kl, nl = int(lin.sum()), int(len(lin)) out["lineages"] = dict(k=kl, n=nl, fraction=kl / nl, expanders=sorted(lin.index[lin == 1]), non_expanders=sorted(lin.index[lin == 0])) # open-ended denominator: how the share moves if capable non-expanders are missing from the list (JUDGEMENT: 10 or 30) out["missing_non_expanders_sensitivity"] = { str(extra): dict(episodes=out["outcome_strict"]["k"] / (out["outcome_strict"]["n"] + extra), lineages=kl / (nl + extra)) for extra in (0, 10, 30)} w = pd.read_csv(os.path.join(CUR, "wave_endings.csv")) out["wave_endings"] = dict( waves=int(len(w)), by_stop=w.stop.value_counts().to_dict(), kept_going_lenient=dict(k=int((w.stop != "retreat").sum()), n=int(len(w)), rule="a stop at a range limit counts as having filled the domain"), kept_going_cautious=dict(k=int((~w.stop.isin(["retreat", "range limit"])).sum()), n=int(len(w)), rule="a stop with empty land beyond the previous longest hop (Lapita) counts as a stall"), every_wave_ended=True) return out, t def illustrative_product(cult, pac, norse, inv, nmc=200000): """Illustrative product of three factors under four explicit scenarios. Conditional on feasibility; valid for the first few relaunch generations in a closed domain; NOT an estimate of the model's expansion probability.""" rng = np.random.default_rng(SEED + 11) ks, ns = cult["outcome_strict"]["k"], cult["outcome_strict"]["n"] kl, nl = cult["outcome_lenient"]["k"], cult["outcome_lenient"]["n"] kL, nL = cult["lineages"]["k"], cult["lineages"]["n"] we = cult["wave_endings"] kp = pac["counts"]["east_persisted"] + (norse["n_frontier_units"] - norse["abandoned"]) np_ = pac["counts"]["east_settled"] + norse["n_frontier_units"] nw = we["waves"] lo_bio = min(inv["stage_summary"]["established_to_spreading"]["gavia_2plus_regions"], 0.25) B = lambda k, n: rng.beta(k + 0.5, n - k + 0.5, nmc) # Beta draws express the sampling noise of the count only q = lambda x: [float(np.percentile(x, p)) for p in (5, 50, 95)] sc = {} # 1. cautious: one faction; lineages with 30 unlisted non-expanders; biological spread rates; settlement-level persistence sc["cautious: one faction; %d of %d lineages plus 30 unlisted non-expanders; biological spread rates 0.10 to %.2f; " "settlement-level persistence" % (kL, nL, lo_bio)] = B(kL, nL + 30) * rng.uniform(0.10, lo_bio, nmc) * B(kp, np_) # 2. lower middle: one faction; lineages as listed; Lapita counted as a stall; mixed persistence c_mix = np.where(rng.random(nmc) < 0.5, B(nw, nw), B(kp, np_)) sc["lower middle: one faction; %d of %d lineages; %d of %d waves kept going; mixed persistence" % (kL, nL, we["kept_going_cautious"]["k"], nw)] = B(kL, nL) * B(we["kept_going_cautious"]["k"], nw) * c_mix # 3. upper middle: three factions; episodes, strict coding; lenient wave coding; mixed persistence a3 = 1 - (1 - B(ks, ns)) ** 3 sc["upper middle: three factions; %d of %d episodes; %d of %d waves kept going; mixed persistence" % (ks, ns, we["kept_going_lenient"]["k"], nw)] = a3 * B(we["kept_going_lenient"]["k"], nw) * c_mix # 4. generous: ten factions; lenient episode coding; lenient wave coding; wave-level persistence a10 = 1 - (1 - B(kl, nl)) ** 10 sc["generous: ten factions; %d of %d episodes; %d of %d waves kept going; wave-level persistence" % (kl, nl, we["kept_going_lenient"]["k"], nw)] = a10 * B(we["kept_going_lenient"]["k"], nw) * B(nw, nw) out = {k: q(v) for k, v in sc.items()} med = [v[1] for v in out.values()] # first-pass style equal mixture, kept for comparison with the first pass (census counts updated) pick = rng.random(nmc) < 0.5 a_one = np.where(pick, B(ks, ns), B(kl, nl)) m = rng.choice([1, 3, 10], nmc) a_civ = 1 - (1 - a_one) ** m b = np.where(rng.random(nmc) < 0.5, B(we["kept_going_lenient"]["k"], nw), rng.uniform(lo_bio, 0.5, nmc)) c = np.where(rng.random(nmc) < 0.5, B(nw, nw), B(kp, np_)) mix = a_civ * b * c return dict(scenarios_q05_q50_q95=out, range_of_scenario_medians=[float(min(med)), float(max(med))], first_pass_style_mixture_q05_q50_q95=q(mix), unit_persistence_counts=dict(persisted=int(kp), n=int(np_)), scope="conditional on feasibility; first few relaunch generations; closed domain of at most about 50 targets", ), np.log10(mix) def floor_sensitivity(log10_mix, nmc=200000): """The first-pass 'prior' multiplied the behavioural product by a feasibility discount drawn log-uniformly between a floor and 1. The floor is a pure JUDGEMENT; this table shows that it sets the answer.""" rng = np.random.default_rng(SEED + 13) lb = rng.choice(log10_mix, nmc) out = {} for k in (2, 4, 6, 10): lp = lb - k * rng.random(nmc) out[f"1e-{k}"] = dict(median=float(10 ** np.median(lp)), q05=float(10 ** np.percentile(lp, 5)), q95=float(10 ** np.percentile(lp, 95)), share_below_1e_3=float((lp < -3).mean()), share_below_1e_10=float((lp < -10).mean())) return out def model_prior_stats(W_range, nmc=400000, seed=5): rng = np.random.default_rng(seed) W = 10 ** rng.uniform(np.log10(W_range[0]), np.log10(W_range[1]), nmc) lp = -W * rng.random(nmc) return dict(W_range=list(W_range), log10_p_quantiles_5_25_50_75_95=[float(np.percentile(lp, p)) for p in (5, 25, 50, 75, 95)], median_p=float(10 ** np.median(lp)), share_below_1e_3=float((lp < -3).mean()), share_above_0p29=float((lp > np.log10(0.29)).mean()), prob_p_gt_1e_3=float((lp > -3).mean()), prob_p_gt_1e_2=float((lp > -2).mean()), prob_p_lt_1e_10=float((lp < -10).mean())), lp def scale_gap(pac, norse): """Orders of magnitude between the analogue and the interstellar case, and between the analogue QUANTITY and the model's parameter.""" hop_med = pac["hops"]["evidence_based_hops_km"]["median"] hop_max = pac["hops"]["evidence_based_hops_km"]["max"] d_star = 4.25 * LY_KM # nearest star system food_py = 2400 * 4184 * 365.25 # J per person-year of food energy # JUDGEMENT: a voyaging canoe and its provisions cost of order 1-10 person-years of labour for a founding party of ~50 polynesian_cost_py_per_settler = (1 / 50, 10 / 50) e_kg = lambda beta: (1 / np.sqrt(1 - beta ** 2) - 1) * C_LIGHT ** 2 tons = (1.0, 100.0) # JUDGEMENT: ship mass per settler, tonnes e_settler = [e_kg(0.01) * 1000 * t for t in tons] # --- the quantity mismatch (revision) --- # Model: p_exp is the chance of becoming a persistent expander that fills a galaxy; model3.py takes the fill time as # max(0.3 Gyr, 1e-4 / v Gyr). JUDGEMENT for the generation count: a front crossing the stellar disc (about 1e5 light-years) # in hops of 5 to 10 light-years makes 1e4 to 2e4 successive relaunch generations. G_lo, G_hi = 1e5 / 10, 1e5 / 5 # Observed 'generation-trials': relaunch steps along the two chains that are documented step by step in this study. gen_pacific = 5 # West Polynesia -> Southern Cooks -> Societies -> Tuamotus -> Mangareva -> Rapa Nui (Ioannidis et al. 2021) gen_norse = 5 # Norway -> Shetland-Orkney -> Faroe -> Iceland -> Greenland -> Vinland n_two = gen_pacific + gen_norse n_gen = 8 * 5 # JUDGEMENT (generous): all eight launched waves credited with five generations each ub = lambda N: float(1 - 0.05 ** (1.0 / N)) # 95% upper bound on a per-trial hazard after N trials with no event need = float(1 - 0.5 ** (1.0 / G_lo)) # hazard that leaves a 50% chance of surviving G_lo generations return dict( pacific_hop_km=dict(median=hop_med, max=hop_max), nearest_star_km=d_star, distance_ratio_log10=[float(np.log10(d_star / hop_max)), float(np.log10(d_star / hop_med))], voyage_time=dict(pacific_weeks=[2, 5], interstellar_years_at_0p01c=4.25 / 0.01, interstellar_years_at_0p1c=4.25 / 0.1, ratio_log10=[float(np.log10(4.25 / 0.1 / (5 / 52))), float(np.log10(4.25 / 0.01 / (2 / 52)))]), energy=dict(kinetic_J_per_kg_at_0p01c=float(e_kg(0.01)), food_J_per_person_year=float(food_py), interstellar_person_years_of_food_energy_per_settler=[float(e / food_py) for e in e_settler], polynesian_person_years_per_settler=list(polynesian_cost_py_per_settler), cost_ratio_log10=[float(np.log10(e_settler[0] / food_py / polynesian_cost_py_per_settler[1])), float(np.log10(e_settler[1] / food_py / polynesian_cost_py_per_settler[0]))]), note="distance gap about 10-11 orders of magnitude; per-settler energy gap about 7-10 orders under the stated assumptions", quantity_mismatch=dict( analogue=dict(targets_in_domain="at most about 50 land units (30 or fewer habitable)", relaunch_generations=[4, 6], wave_duration_yr=[250, 800], outcome="every documented wave ended"), model_parameter=dict(definition="persistent expander that fills a galaxy (model3.py: fill time at least 0.3 Gyr)", targets="about 1e11 stars", relaunch_generations_along_front=[G_lo, G_hi], duration_yr=3e8), log10_gap=dict(targets=float(np.log10(1e11 / 50)), generations=[float(np.log10(G_lo / 6)), float(np.log10(G_hi / 4))], duration=[float(np.log10(3e8 / 800)), float(np.log10(3e8 / 250))])), hazard_bound=dict( generation_trials_two_documented_chains=n_two, generation_trials_generous=n_gen, upper95_per_generation_stopping_hazard_if_no_stop_was_intrinsic=dict(two_chains=ub(n_two), generous=ub(n_gen)), hazard_allowing_50pct_survival_over_1e4_generations=need, generation_trials_needed_to_bound_hazard_that_low=float(np.log(0.05) / np.log(1 - need)), note="even if none of the observed stops reflected an intrinsic hazard, at most about 40 observed wave-generations bound the " "per-generation hazard only below about 0.07; sustaining 1e4 generations needs it below about 7e-5. The analogues are " "therefore silent, in both directions, about the long-run factor that dominates the model's parameter."), ) def run_model(verbose=True): """Re-run the project's integrated model (frozen copy) with several ASSUMED expansion priors. A sensitivity test: it shows what the model would say IF the prior were changed, not that the prior should be changed.""" vend = os.path.join(HERE, "vendor_model3") if not os.path.exists(os.path.join(vend, "model3.py")): return dict(error="vendor_model3 not present") sys.path.insert(0, vend) sys.dont_write_bytecode = True import model3 as M runs = {"base W=(1,30)": {}, "log-uniform 1e-2 to 1, W=(2,2)": dict(W_exp=(2.0, 2.0)), "log-uniform 1e-4 to 1, W=(4,4)": dict(W_exp=(4.0, 4.0)), "log-uniform 1e-6 to 1, W=(6,6)": dict(W_exp=(6.0, 6.0)), "first-pass branch W=(1,3)": dict(W_exp=(1.0, 3.0)), "first-pass branch W=(2,5)": dict(W_exp=(2.0, 5.0))} out = {} for name, kw in runs.items(): r = M.run(M.vary(M.BASE, n=200000, **kw)) out[name] = dict(prior=r["prior"], schools={s: dict(posterior=r["schools"][s]["posterior"], evidence=r["schools"][s]["evidence"], ess=r["schools"][s]["ess"], p_no_other_ever_in_galaxy=r["schools"][s]["p_no_other_ever_in_galaxy"]) for s in ("site", "ssa", "sia")}) base = out["base W=(1,30)"]["schools"] for name in list(out): if name.startswith("base"): continue for s in ("site", "ssa"): out[name]["schools"][s]["bayes_factor_base_over_this"] = base[s]["evidence"] / out[name]["schools"][s]["evidence"] # equal mixture of the base prior and the first-pass W=(2,5) branch, kept from the first pass for comparison alt = out["first-pass branch W=(2,5)"]["schools"] mix = {} for s in ("site", "ssa"): zb, za = base[s]["evidence"], alt[s]["evidence"] wb, wa = 0.5 * zb, 0.5 * za mix[s] = dict(posterior={c: (wb * base[s]["posterior"][c] + wa * alt[s]["posterior"][c]) / (wb + wa) for c in M.CLASSES}) out["equal mixture of base and W=(2,5)"] = dict(schools=mix) out["class_labels"] = M.CLASS_LABELS if verbose: for k, v in out.items(): if k == "class_labels": continue p = v["schools"]["site"]["posterior"] print(f" {k:36s} site: " + " ".join(f"{c}={p[c]:.3f}" for c in M.CLASSES[:6])) return out def run(pac, norse, inv, verbose=True, with_model=True): cult, tab = cultures() ill, lmix = illustrative_product(cult, pac, norse, inv) priors = {} lp = {} for name, W in (("base W=(1,30)", (1.0, 30.0)), ("W=(2,2)", (2.0, 2.0)), ("W=(4,4)", (4.0, 4.0)), ("W=(6,6)", (6.0, 6.0)), ("W=(1,3)", (1.0, 3.0)), ("W=(2,5)", (2.0, 5.0))): priors[name], lp[name] = model_prior_stats(W) floors = floor_sensitivity(lmix) res = dict(cultures=cult, illustrative_product=ill, feasibility_floor_sensitivity=floors, model_priors=priors, scale_gap=scale_gap(pac, norse)) res["what_can_be_said"] = dict( conditional_short_run=("given feasibility and a capable culture, launching a wave and carrying it through the first few " "relaunch generations of a small closed domain was not rare: scenario medians %.2f to %.2f" % tuple(ill["range_of_scenario_medians"])), long_run="the analogues cannot constrain the probability of sustaining expansion beyond a few generations, in either direction", model="no change to the project model's W_exp is recommended on the strength of this study; the model re-runs are a " "sensitivity test of assumed priors") if with_model: res["model_rerun"] = run_model(verbose) rng = np.random.default_rng(SEED + 14) lbf = {k: (rng.choice(lmix, 50000) - int(k[3:]) * rng.random(50000)) for k in ("1e-2", "1e-4", "1e-6")} np.savez_compressed(os.path.join(DATA, "prior_samples.npz"), log10_product_mixture=lmix[:50000], log10_p_base=lp["base W=(1,30)"][:50000], log10_floor_1e_2=lbf["1e-2"], log10_floor_1e_4=lbf["1e-4"], log10_floor_1e_6=lbf["1e-6"]) dump_json(res, "synthesis_results.json") if verbose: print("Maritime census: episodes strict %d/%d, lenient %d/%d; lineages %d/%d" % ( cult["outcome_strict"]["k"], cult["outcome_strict"]["n"], cult["outcome_lenient"]["k"], cult["outcome_lenient"]["n"], cult["lineages"]["k"], cult["lineages"]["n"])) for k, v in ill["scenarios_q05_q50_q95"].items(): print(" %-60s %s" % (k[:60], [round(x, 3) for x in v])) print(" floor sensitivity:", {k: (round(v["median"], 4), round(v["share_below_1e_3"], 2)) for k, v in floors.items()}) return res if __name__ == "__main__": import json pac = json.load(open(os.path.join(DATA, "pacific_results.json"))) norse = json.load(open(os.path.join(DATA, "norse_results.json"))) inv = json.load(open(os.path.join(DATA, "invasion_results.json"))) run(pac, norse, inv)