"""Pacific settlement wave: hop distances, timing, abandonment, and a percolation fit on the real island graph. Revised after the referee report. Changes that matter: * reach: the share of targets settled is reported under several denominators, with the outcome-dependent one flagged; * abandonment: an island-level analysis (persisted archipelagos split into islands) with unit-clustered inference and an aggregation-radius profile sits beside the unit-level regression; * percolation: a second node set that does not depend on the outcome (all land units, with a fitted habitability term); * timing: 'R' is called what it is, a growth factor of the count of settled units, with subset sensitivity.""" import heapq import json import os import numpy as np import pandas as pd from common import CUR, DATA, SEED, dump_json, haversine, jeffreys, logistic_fit, logistic_fit_fast GENERIC_WINDOWS = {(1190, 1293), (1300, 1600)} # short-chronology windows that are not specific to the unit W11_SPECIFIC = ["Society Islands", "Marquesas", "Mangareva (Gambier)", "Rapa Nui", "Hawaiian Islands", "New Zealand", "Southern Cooks"] CORE_SOURCES = ["Fiji", "Tonga", "Samoa", "Uvea", "Futuna"] # Lapita-period units; the other West units have no firm date DOUBTFUL = ["Howland", "Necker", "Auckland Islands"] # abandoned units whose evidence speaks against a resident settlement SOUTH_AMERICA_COAST = (-27.1, -70.9) # nearest stretch of the Chilean coast to Rapa Nui (Atacama coast) def load(): geo = json.load(open(os.path.join(DATA, "wikidata_geo_cache.json"))) a = pd.read_csv(os.path.join(CUR, "pacific_anchor_islands.csv")) u = pd.read_csv(os.path.join(CUR, "pacific_units.csv")) lat, lon = [], [] for _, r in a.iterrows(): if isinstance(r.geo_key, str): c = geo[r.geo_key]["wd_coord"] lat.append(c[0]); lon.append(c[1]) else: lat.append(r.manual_lat); lon.append(r.manual_lon) a["lat"], a["lon"] = lat, lon assert set(a.unit) == set(u.unit) units = list(u.unit) n = len(units) D = np.zeros((n, n)) for i, ui in enumerate(units): ai = a[a.unit == ui] for j in range(i + 1, n): aj = a[a.unit == units[j]] D[i, j] = D[j, i] = haversine(ai.lat.values[:, None], ai.lon.values[:, None], aj.lat.values[None, :], aj.lon.values[None, :]).min() cen = a.groupby("unit")[["lat", "lon"]].first().reindex(units) # first anchor as plotting position u["lat"], u["lon"] = cen.lat.values, cen.lon.values u["settled"] = u.status.isin(["persisted", "abandoned"]) return a, u, D def bottleneck_from_sources(D, node_ok, sources): """Minimax-path ('bottleneck') distance from a set of source nodes to every allowed node: the smallest single-hop range that connects the node to the sources through allowed nodes.""" n = len(D) best = np.full(n, np.inf) prev = np.full(n, -1) pq = [] for s in sources: best[s] = 0.0 heapq.heappush(pq, (0.0, s)) while pq: b, i = heapq.heappop(pq) if b > best[i]: continue for j in range(n): if j == i or not node_ok[j]: continue nb = max(b, D[i, j]) if nb < best[j]: best[j] = nb prev[j] = i heapq.heappush(pq, (nb, j)) return best, prev def single_linkage(D, idx, L): """Connected components of the graph on nodes idx with edges shorter than L.""" idx = list(idx) lab = {i: k for k, i in enumerate(idx)} changed = True while changed: changed = False for a_ in idx: for b_ in idx: if a_ < b_ and D[a_, b_] < L and lab[a_] != lab[b_]: old, new = lab[b_], lab[a_] for k in idx: if lab[k] == old: lab[k] = new changed = True return lab def sample_dates(u, chron, rng, nmc): lo = u[f"{chron}_lo"].values.astype(float) hi = u[f"{chron}_hi"].values.astype(float) return lo[None, :] + (hi - lo)[None, :] * rng.random((nmc, len(u))) def timing(u, D, chron, rng, nmc=20000, subset=None): """Monte-Carlo over date windows: pause, pulse duration, growth rate of the count of settled units, relaunch delays. subset: optional boolean mask over units restricting which dated East units enter (sensitivity runs).""" t = sample_dates(u, chron, rng, nmc) names = list(u.unit) ix = {k: i for i, k in enumerate(names)} east = (u.region == "E").values & u.settled.values & np.isfinite(t[0]) if subset is not None: east = east & subset te = t[:, east] first = te.min(1) pause = first - t[:, ix["Samoa"]] srt = np.sort(te, axis=1) ne = te.shape[1] k10, k50, k90 = int(np.ceil(0.1 * ne)) - 1, int(np.ceil(0.5 * ne)) - 1, int(np.ceil(0.9 * ne)) - 1 t10, t50, t90 = srt[:, k10], srt[:, k50], srt[:, k90] # exponential-phase growth rate: from the first unit to the time half of the dated units exist r_half = np.log(k50 + 1) / np.maximum(t50 - first, 1.0) # burst-only rate: from the second unit to the time 90% of the dated units exist (drops the solitary first unit) r_burst = np.log((k90 + 1) / 2.0) / np.maximum(t90 - srt[:, 1], 1.0) # relaunch delays for evidence-based parent -> child pairs lags = {} for i, r in u.iterrows(): if isinstance(r.parent, str) and r.region == "E" and np.isfinite(t[0, i]) and np.isfinite(t[0, ix[r.parent]]): lags[(r.parent, r.unit)] = t[:, i] - t[:, ix[r.parent]] parents = sorted({p for p, _ in lags}) relaunch = {p: np.min(np.stack([v for (pp, c), v in lags.items() if pp == p]), axis=0) for p in parents} rel_mat = np.stack([relaunch[p] for p in parents], axis=1) east_parents = [k for k, p in enumerate(parents) if u.region[ix[p]] == "E"] rel_east = rel_mat[:, east_parents] grid = np.arange(-1200, 1801, 10) cum = (te[:, :, None] <= grid[None, None, :]).sum(1) q = lambda x: [float(np.percentile(x, p)) for p in (2.5, 16, 50, 84, 97.5)] out = dict( chronology=chron, n_dated_east_units=int(ne), first_east_settlement=q(first), pause_years=q(pause), t10=q(t10), t50=q(t50), t90=q(t90), pulse_10_90_years=q(t90 - t10), growth_rate_first_to_half_per_yr=q(r_half), doubling_time_first_to_half_yr=q(np.log(2) / r_half), growth_rate_burst_per_yr=q(r_burst), doubling_time_burst_yr=q(np.log(2) / r_burst), R_per_30yr_generation_first_to_half=q(np.exp(30 * r_half)), R_per_30yr_generation_burst=q(np.exp(30 * r_burst)), new_units_after_1350=q((te > 1350).sum(1)), units_by_1350=q((te <= 1350).sum(1)), pair_lag_median={f"{p} -> {c}": float(np.median(v)) for (p, c), v in lags.items()}, pair_lag_p16_p84={f"{p} -> {c}": [float(np.percentile(v, 16)), float(np.percentile(v, 84))] for (p, c), v in lags.items()}, pair_lag_frac_negative={f"{p} -> {c}": float((v < 0).mean()) for (p, c), v in lags.items()}, relaunch_delay_by_parent_median={p: float(np.median(relaunch[p])) for p in parents}, relaunch_delay_east_parents_median_of_medians=float(np.median([np.median(relaunch[parents[k]]) for k in east_parents])), relaunch_delay_east_parents_pooled=q(rel_east.ravel()), relaunch_delay_east_parents_frac_under_100yr=float((rel_east < 100).mean()), ) curve = pd.DataFrame(dict(year=grid, lo=np.percentile(cum, 16, axis=0), med=np.percentile(cum, 50, axis=0), hi=np.percentile(cum, 84, axis=0))) return out, curve def _sigmoid(z): return 1.0 / (1.0 + np.exp(-z)) def percolation_mc(D, settled0, origin, must_settle, must_not, p_grid, s_grid, R, rng, nsim=4000, extra=None): """Landis-type site percolation on the real island graph. A settled node is a 'coloniser' with probability p (the origin always is). A coloniser reaches each unsettled node within range R independently with probability s (s = 1 is Landis's original rule). Returns, on the (p, s) grid, the probability of cover: settling every node in must_settle (the must_not nodes are not allowed to act as stepping stones); exact: that, and leaving every must_not node empty; frac : the mean fraction of must_settle nodes reached (a direct, non-asymptotic measure of how far the wave gets). 'exact' is computed without extra Monte-Carlo noise: given the final set of active colonisers, a must_not node with k of them in range stays empty with probability (1 - s)**k (its incoming links are independent of the rest). extra (optional) switches on the outcome-independent variant ("model B"): dict with D_extra : distances from the simulation nodes to further land units that were never settled, x_set, z_set : log10 area and subantarctic-or-rock flag of the must_settle nodes, x_un, z_un : the same for the unsettled units (the must_not nodes first, then the D_extra columns). Every land unit is then habitable with probability h = logistic(a + b x + c z); an unsettled unit is either not habitable or habitable and missed, so the likelihood of the observed pattern is prod_settled h_i x E[ 1{all settled units reached} prod_unsettled (1 - h_j + h_j (1 - s)**k_j) ], maximised over (a, b, c) in every grid cell (Nelder-Mead; the expectation uses at most 800 covering simulations).""" from scipy.optimize import minimize n = len(D) adj = (D <= R) & ~np.eye(n, dtype=bool) adj_dyn = adj.copy() adj_dyn[:, must_not] = False cover = np.zeros((len(p_grid), len(s_grid))) exact = np.zeros_like(cover) frac = np.zeros_like(cover) likB = np.zeros_like(cover) parB = np.full(cover.shape + (3,), np.nan) if extra is not None: adjU = np.column_stack([adj[:, must_not], extra["D_extra"] <= R]).astype(np.int32) xs, zs, xu, zu = extra["x_set"], extra["z_set"], extra["x_un"], extra["z_un"] for a_, p in enumerate(p_grid): for b_, s in enumerate(s_grid): col = rng.random((nsim, n)) < p col[:, origin] = True edge = (rng.random((nsim, n, n)) < s) & adj_dyn[None] settled = np.tile(settled0, (nsim, 1)) for _ in range(n): active = settled & col new = (active[:, :, None] & edge).any(1) & ~settled if not new.any(): break settled |= new ok = settled[:, must_settle].all(1) active = settled & col k = (active[:, :, None] & adj[:, must_not][None]).sum(1) # active colonisers in range of each must_not node w = np.prod(np.power(1.0 - s, k), axis=1) # 0**0 = 1: nobody in range, stays empty cover[a_, b_] = ok.mean() exact[a_, b_] = (ok * w).mean() frac[a_, b_] = settled[:, must_settle].mean() if extra is not None and ok.any(): sel = np.where(ok)[0][:800] K = active[sel].astype(np.int32) @ adjU # (sims, unsettled units) q = np.power(1.0 - s, K) # chance each unsettled unit is missed if habitable pok = ok.mean() def nll(th): hs = _sigmoid(th[0] + th[1] * xs + th[2] * zs) hu = _sigmoid(th[0] + th[1] * xu + th[2] * zu) m = np.mean(np.prod(1.0 - hu[None, :] + hu[None, :] * q, axis=1)) return -(np.sum(np.log(np.maximum(hs, 1e-300))) + np.log(max(pok * m, 1e-300))) best = minimize(nll, extra.get("theta0", (-0.5, 1.3, -3.0)), method="Nelder-Mead", options=dict(xatol=0.02, fatol=0.002, maxiter=400)) likB[a_, b_] = float(np.exp(-best.fun)) parB[a_, b_] = best.x return cover, exact, frac, likB, parB def island_level(u, a, D_units, verbose=False): """Abandonment at island grain (revision). The unit table merges the persisted small islands into archipelagos, while every abandoned island is its own unit; here persisted archipelagos are split into islands. Inference is clustered on the 29 settlement units, because islands of one archipelago share their fate.""" from scipy.stats import chi2 import island_sources isl = pd.read_csv(os.path.join(CUR, "pacific_islands.csv"), keep_default_na=False) cache = island_sources.load_cache() isl["lat"] = [cache[t]["lat"] for t in isl.wp_title] isl["lon"] = [cache[t]["lon"] for t in isl.wp_title] isl["land_km2"] = isl.land_km2.astype(float) # every other piece of land in the unit table (West Polynesia, never-settled units, eastern rim) as extra land points covered = set(isl.unit) geo = json.load(open(os.path.join(DATA, "wikidata_geo_cache.json"))) ext = [] for unit, grp in a[~a.unit.isin(covered)].groupby("unit"): tot = float(u.land_km2[u.unit == unit].iloc[0]) claims = [] for _, r in grp.iterrows(): c = geo.get(r.geo_key, {}).get("area_claims") if isinstance(r.geo_key, str) else None claims.append(float(c[0][0]) if c else np.nan) claims = np.array(claims) # cap the Wikidata claims at the unit total (several atoll claims include the lagoon); share the rest equally if np.isfinite(claims).all() and claims.sum() <= 1.25 * tot: areas = claims else: areas = np.full(len(grp), tot / len(grp)) for (_, r), ar in zip(grp.iterrows(), areas): ext.append(dict(island=r.island, unit=unit, lat=r.lat, lon=r.lon, land_km2=float(ar))) ext = pd.DataFrame(ext) lat = np.concatenate([isl.lat.values, ext.lat.values]).astype(float) lon = np.concatenate([isl.lon.values, ext.lon.values]).astype(float) land = np.concatenate([isl.land_km2.values, ext.land_km2.values]).astype(float) n = len(isl) Dall = haversine(lat[:n, None], lon[:n, None], lat[None, :], lon[None, :]) Dall[np.arange(n), np.arange(n)] = np.inf Dss = Dall[:, :n] isl["nn_settled_km"] = Dss.min(1) # nearest other island with settlement evidence isl["n_settled_300km"] = (Dss < 300).sum(1) isl["iso_big_km"] = np.where(land[None, :] >= 100, Dall, np.inf).min(1) # nearest other land of at least 100 km2 (any status) isl["land_within_300km"] = (np.where(Dall < 300, land[None, :], 0.0)).sum(1) radii = [0, 50, 100, 200, 300, 500, 1000] for r_ in radii: isl[f"cluster_land_{r_}km"] = isl.land_km2.values + (np.where(Dall < r_, land[None, :], 0.0)).sum(1) isl.to_csv(os.path.join(DATA, "pacific_islands_derived.csv"), index=False) y_all = (isl.contact == "abandoned").values.astype(float) doubt_p = (isl.contact == "doubtful").values doubt_a = isl.unit.isin(DOUBTFUL).values & (y_all == 1) samples = { "main (doubtful persisted islands dropped)": (~doubt_p, y_all), "doubtful persisted islands kept as persisted": (np.ones(n, bool), y_all), "doubtful persisted islands counted as abandoned": (np.ones(n, bool), np.where(doubt_p, 1.0, y_all)), "main, three doubtful abandoned settlements dropped": (~doubt_p & ~doubt_a, y_all), } covs = { "iso_big_km": lambda d: np.log10(d.iso_big_km.values), "nn_settled_km": lambda d: np.log10(d.nn_settled_km.values), "n_settled_300km": lambda d: np.log10(1.0 + d.n_settled_300km.values), "land_within_300km": lambda d: np.log10(1.0 + d.land_within_300km.values), } rng = np.random.default_rng(SEED + 21) def one_fit(X, y, ridge=0.0): b, cov, ll = logistic_fit(X, y, ridge=ridge) return b, np.sqrt(np.diag(cov)), ll def fit_block(mask, y, cov_name): d = isl[mask].reset_index(drop=True) yy = y[mask] X = np.column_stack([np.log10(d.land_km2.values), covs[cov_name](d)]) b, se, ll = one_fit(X, yy) _, _, ll_a = one_fit(X[:, :1], yy) _, _, ll_i = one_fit(X[:, 1:], yy) # cluster bootstrap: resample settlement units with replacement units = d.unit.values ulist = sorted(set(units)) idx = {k: np.where(units == k)[0] for k in ulist} bs = [] for _ in range(1000): pick = rng.integers(0, len(ulist), len(ulist)) j = np.concatenate([idx[ulist[k]] for k in pick]) if yy[j].min() == yy[j].max(): continue bs.append(logistic_fit_fast(X[j], yy[j], ridge=0.05)) bs = np.array(bs) # one island per unit: removes pseudo-replication altogether (n = number of units) one = [] for _ in range(400): j = np.array([rng.choice(idx[k]) for k in ulist]) one.append(logistic_fit_fast(X[j], yy[j], ridge=0.05)) one = np.array(one) return dict( n=int(len(yy)), n_abandoned=int(yy.sum()), n_units=int(len(ulist)), coef=dict(log10_area=float(b[1]), isolation=float(b[2])), naive_se=dict(log10_area=float(se[1]), isolation=float(se[2])), naive_lr_p=dict(log10_area=float(chi2.sf(2 * (ll - ll_i), 1)), isolation=float(chi2.sf(2 * (ll - ll_a), 1))), cluster_boot95=dict(log10_area=[float(np.percentile(bs[:, 1], 2.5)), float(np.percentile(bs[:, 1], 97.5))], isolation=[float(np.percentile(bs[:, 2], 2.5)), float(np.percentile(bs[:, 2], 97.5))]), cluster_boot_frac_wrong_sign=dict(log10_area=float((bs[:, 1] > 0).mean()), isolation=float((np.sign(bs[:, 2]) != np.sign(b[2])).mean())), one_island_per_unit=dict(median_coef=dict(log10_area=float(np.median(one[:, 1])), isolation=float(np.median(one[:, 2]))), p2p5_p97p5=dict(log10_area=[float(np.percentile(one[:, 1], 2.5)), float(np.percentile(one[:, 1], 97.5))], isolation=[float(np.percentile(one[:, 2], 2.5)), float(np.percentile(one[:, 2], 97.5))])), ) out = dict(n_islands=int(n), n_abandoned=int(y_all.sum()), n_inhabited=int((isl.contact == "inhabited").sum()), n_doubtful=int(doubt_p.sum()), completeness_note="Tuamotus represented by 15 of 76 atolls; other persisted archipelagos by every island with an infobox land area", area_check_entered_vs_infobox=island_sources.check_areas()) m_main, y_main = samples["main (doubtful persisted islands dropped)"] d = isl[m_main] small = d.land_km2.values < 100 out["share_abandoned_under_100km2"] = dict(abandoned=int(y_main[m_main][small].sum()), islands=int(small.sum()), fraction=float(y_main[m_main][small].mean())) cls = [(0, 10), (10, 100), (100, 1000), (1000, 1e9)] out["by_area_class"] = [dict(area_km2=f"{lo}-{hi if hi < 1e8 else 'inf'}", islands=int(((d.land_km2 >= lo) & (d.land_km2 < hi)).sum()), abandoned=int((y_main[m_main] * ((d.land_km2 >= lo) & (d.land_km2 < hi))).sum())) for lo, hi in cls] sm = d[small] ysm = y_main[m_main][small] from scipy.stats import mannwhitneyu out["small_islands_medians"] = { c: dict(abandoned=float(np.median(sm[c][ysm == 1])), persisted=float(np.median(sm[c][ysm == 0])), mannwhitney_p_naive=float(mannwhitneyu(sm[c][ysm == 1], sm[c][ysm == 0]).pvalue)) for c in ("land_km2", "nn_settled_km", "iso_big_km", "n_settled_300km", "land_within_300km")} out["fits"] = {sn: {cn: fit_block(m, y, cn) for cn in covs} for sn, (m, y) in samples.items()} # aggregation-radius profile: which grain of 'area' predicts abandonment best? prof = [] units = isl.unit.values[m_main] ulist = sorted(set(units)) yy = y_main[m_main] for r_ in radii: x = np.log10(isl[f"cluster_land_{r_}km"].values[m_main])[:, None] b, se, ll = one_fit(x, yy) _, _, ll0 = one_fit(np.zeros((len(yy), 0)), yy) # leave-one-unit-out prediction pl = np.zeros(len(yy)) for k in ulist: te = units == k bb = logistic_fit_fast(x[~te], yy[~te], ridge=0.05) pl[te] = _sigmoid(bb[0] + x[te, 0] * bb[1]) pos, neg = pl[yy == 1], pl[yy == 0] auc = float(((pos[:, None] > neg[None, :]).sum() + 0.5 * (pos[:, None] == neg[None, :]).sum()) / (len(pos) * len(neg))) bs = [] idx = {k: np.where(units == k)[0] for k in ulist} for _ in range(500): pick = rng.integers(0, len(ulist), len(ulist)) j = np.concatenate([idx[ulist[k]] for k in pick]) if yy[j].min() == yy[j].max(): continue bs.append(logistic_fit_fast(x[j], yy[j], ridge=0.05)[1]) prof.append(dict(radius_km=r_, slope_per_dex=float(b[1]), naive_se=float(se[1]), cluster_boot95=[float(np.percentile(bs, 2.5)), float(np.percentile(bs, 97.5))], loglik=float(ll), loglik_null=float(ll0), deviance_explained=float(1 - ll / ll0), leave_one_unit_out_auc=auc, land_at_p50_km2=float(10 ** (-b[0] / b[1])))) out["aggregation_radius_profile"] = prof best = max(prof, key=lambda r: r["loglik"]) out["aggregation_radius_best"] = dict(radius_km=best["radius_km"], loglik_gain_over_island_area=float(best["loglik"] - prof[0]["loglik"])) if verbose: print(json.dumps(out, indent=1, default=float)[:5000]) return out, isl def run(verbose=True): rng = np.random.default_rng(SEED) a, u, D = load() names = list(u.unit) ix = {k: i for i, k in enumerate(names)} n = len(u) res = {} settled = u.settled.values W = (u.region == "W").values E = (u.region == "E").values Es = E & settled doubtful = u.unit.isin(DOUBTFUL).values # ------------------------------------------------------------------ counts res["counts"] = dict( n_units=int(n), n_anchor_islands=int(len(a)), by_region_status=u.groupby(["region", "status"]).size().rename("n").reset_index().to_dict("records"), east_settled=int(Es.sum()), east_persisted=int((E & (u.status == "persisted")).sum()), east_abandoned=int((E & (u.status == "abandoned")).sum()), east_traces=int((E & (u.status == "traces")).sum()), east_none=int((E & (u.status == "none")).sum()), verification_grades=u.grade.value_counts().to_dict(), abandoned_evidence=u.evidence[E & (u.status == "abandoned").values].value_counts().to_dict(), doubtful_settlements=DOUBTFUL, ) # ------------------------------------------------------------------ hop distances Dm = D.copy() np.fill_diagonal(Dm, np.inf) nearest_settled = np.where(settled[None, :], Dm, np.inf).min(1) nearest_settled_name = [names[j] for j in np.where(settled[None, :], Dm, np.inf).argmin(1)] big = (u.land_km2.values >= 100) iso_big = np.where(big[None, :], Dm, np.inf).min(1) pers = (u.status == "persisted").values iso_pers = np.where(pers[None, :], Dm, np.inf).min(1) sources = [i for i in range(n) if W[i]] bn_all, prev_all = bottleneck_from_sources(D, settled, sources) bn_pers, _ = bottleneck_from_sources(D, pers, sources) u["nearest_settled_km"] = nearest_settled u["nearest_settled_unit"] = nearest_settled_name u["isolation_big_km"] = iso_big u["isolation_persisted_km"] = iso_pers u["bottleneck_km"] = np.where(settled, bn_all, np.nan) u["bottleneck_without_abandoned_km"] = np.where(pers, bn_pers, np.nan) u["bottleneck_via"] = [names[p] if p >= 0 else "" for p in prev_all] u["parent_km"] = [D[i, ix[r.parent]] if isinstance(r.parent, str) else np.nan for i, r in u.iterrows()] # potential parents within a few candidate ranges for R in (2100, 2600, 3500): u[f"n_settled_within_{R}"] = [(settled & (Dm[i] <= R)).sum() for i in range(n)] hop_min = nearest_settled[Es] ev = u[E & settled & u.parent.notna()] sa = float(haversine(u.lat[ix["Rapa Nui"]], u.lon[ix["Rapa Nui"]], *SOUTH_AMERICA_COAST)) D_net = float(np.nanmax(u.bottleneck_km[Es])) D_net_unit = u.unit[Es].values[int(np.nanargmax(u.bottleneck_km[Es].values))] corner = {k: dict(bottleneck_km=float(u.bottleneck_km[ix[k]]), via=u.bottleneck_via[ix[k]], bottleneck_without_abandoned_km=float(u.bottleneck_without_abandoned_km[ix[k]])) for k in ("Hawaiian Islands", "Rapa Nui", "New Zealand")} longest_ev = ev.sort_values("parent_km").iloc[-1] # sensitivity: sources restricted to the dated Lapita-period core; the undated West units join the pulse domain core = np.array([k in CORE_SOURCES for k in names]) bn_core, _ = bottleneck_from_sources(D, settled, [i for i in range(n) if core[i]]) dom_core = settled & ~core # sensitivity: the three doubtful settlements not available as stepping stones bn_nd, _ = bottleneck_from_sources(D, settled & ~doubtful, sources) res["hops"] = dict( definition="distances between units are the minimum great-circle distance between their anchor islands", nearest_settled_neighbour_km=dict(n=int(Es.sum()), median=float(np.median(hop_min)), q25=float(np.percentile(hop_min, 25)), q75=float(np.percentile(hop_min, 75)), max=float(hop_min.max()), max_unit=u.unit[Es].values[int(hop_min.argmax())]), evidence_based_hops_km=dict(n=int(len(ev)), median=float(ev.parent_km.median()), q25=float(ev.parent_km.quantile(.25)), q75=float(ev.parent_km.quantile(.75)), max=float(ev.parent_km.max()), max_pair=f"{longest_ev.parent} -> {longest_ev.unit}", pairs={f"{r.parent} -> {r.unit}": round(float(r.parent_km)) for _, r in ev.iterrows()}), network_bottleneck_km=D_net, network_bottleneck_unit=D_net_unit, network_bottleneck_note="defined from the settled set: the smallest single-hop range that connects every settled unit to " "West Polynesia through settled units; it is a property of the outcome, not a measured capability", corners=corner, rapa_nui_to_south_america_km=sa, attested_range_summary=dict( network_needs_km=[D_net, float(max(c["bottleneck_without_abandoned_km"] for c in corner.values()))], longest_documented_settlement_hop_km=float(ev.parent_km.max()), south_america_contact_km=sa, unsettled_later_habitable_beyond_network={k: float(nearest_settled[ix[k]]) for k in ("Juan Fernandez", "Galapagos")}), core_sources_sensitivity=dict(sources=CORE_SOURCES, network_bottleneck_km=float(np.nanmax(bn_core[dom_core])), bottleneck_unit=names[int(np.where(dom_core, bn_core, -1).argmax())], domain_units=int(dom_core.sum())), without_doubtful_settlements=dict(network_bottleneck_km=float(np.nanmax(bn_nd[Es & ~doubtful]))), frac_evidence_hops_below_half_bottleneck=float((ev.parent_km < 0.5 * D_net).mean()), frac_nearest_hops_below_half_bottleneck=float((hop_min < 0.5 * D_net).mean()), hop_over_range_quantiles={"range_used_km": D_net, "nearest_neighbour": [float(np.percentile(hop_min / D_net, p)) for p in (25, 50, 75, 100)], "evidence_based": [float(ev.parent_km.quantile(p) / D_net) for p in (.25, .5, .75, 1.0)]}, ) # ------------------------------------------------------------------ timing res["timing"] = {} curves = [] for chron in ("S", "A"): out, curve = timing(u, D, chron, rng) res["timing"][chron] = out curve["chronology"] = chron curves.append(curve) pd.concat(curves).to_csv(os.path.join(DATA, "pacific_cumulative_settlement_curves.csv"), index=False) # relaunch delays along the genetic settlement chain (point dates, Ioannidis et al. 2021) gc = pd.read_csv(os.path.join(CUR, "genetic_settlement_chain_ioannidis2021.csv")) gc["interval_yr"] = gc.child_date_CE - gc.parent_date_CE first_daughter = gc.groupby("parent_island").interval_yr.min() # delay until a settlement founds its first daughter res["timing"]["genetic_chain"] = dict( source="Ioannidis et al. 2021, point dates (no intervals propagated)", links={f"{r.parent_island} -> {r.child_island}": int(r.interval_yr) for _, r in gc.iterrows()}, relaunch_delay_first_daughter_yr={k: int(v) for k, v in first_daughter.items()}, relaunch_delay_median_yr=float(first_daughter.median()), relaunch_delay_median_excluding_first_step_yr=float(first_daughter.drop("Rarotonga").median()), chain_rarotonga_to_rapa_nui=dict(years=int(1210 - 830), steps=4, mean_years_per_step=(1210 - 830) / 4), caveat="four point dates from the genetic (alternative) chronology; under the short radiocarbon chronology the Southern " "Cooks (1231-1290) post-date their supposed daughter the Society Islands (1025-1121), so the delay cannot be " "combined with short-chronology dates without mixing chronologies", ) # growth factor per other intervals: G(tau) = exp(r * tau) = 2 ** (tau / doubling time) tau_gen = res["timing"]["genetic_chain"]["relaunch_delay_median_yr"] for chron in ("S", "A"): td = np.array(res["timing"][chron]["doubling_time_first_to_half_yr"]) res["timing"][chron]["R_per_relaunch_delay_first_to_half"] = dict( tau_yr=tau_gen, p2p5_p16_p50_p84_p97p5=[float(x) for x in 2 ** (tau_gen / td[::-1])], mixes_chronologies=bool(chron == "S")) # late phase: new units founded after 1350 per unit existing in 1350 (per remaining ~420 years to 1770) for chron in ("S", "A"): t_ = res["timing"][chron] t_["late_new_units_per_existing_unit"] = float(t_["new_units_after_1350"][2] / max(t_["units_by_1350"][2], 1)) # sensitivity of the growth summary to WHICH units carry dates (separate random stream; the main results are unchanged) rng_t = np.random.default_rng(SEED + 31) generic = np.array([(r.S_lo, r.S_hi) in GENERIC_WINDOWS for _, r in u.iterrows()]) subsets = {"all dated units": None, "units with a unit-specific window": ~generic, "the seven archipelago windows of Wilmshurst et al.": u.unit.isin(W11_SPECIFIC).values} sens = {} for label, msk in subsets.items(): for chron in ("S", "A"): o, _ = timing(u, D, chron, rng_t, nmc=8000, subset=msk) td = o["doubling_time_first_to_half_yr"] sens[f"{label} | {chron}"] = dict( n_dated=o["n_dated_east_units"], doubling_time_yr=[td[0], td[2], td[4]], growth_factor_per_30yr=[o["R_per_30yr_generation_first_to_half"][k] for k in (0, 2, 4)], growth_factor_per_relaunch_delay=[float(2 ** (tau_gen / td[k])) for k in (4, 2, 0)], pulse_10_90_years=[o["pulse_10_90_years"][k] for k in (0, 2, 4)]) res["timing"]["growth_sensitivity"] = sens res["timing"]["generic_window_units"] = list(u.unit[generic & Es]) res["timing"]["growth_note"] = ("growth factor = exp(r x interval) with r = ln(number of dated units at the half-way point) / " "(time from the first unit to the half-way point); it restates the doubling time of the COUNT of " "dated units and is not a measured number of daughters per settlement. Under the short chronology " "the curve is a step (one archipelago, then one pulse), not an exponential.") # documented family tree: how many settlements are known to have founded others? par = u[E & settled & u.parent.notna()].parent.value_counts() par_e = par[[k for k in par.index if u.region[ix[k]] == "E"]] n_e = int(Es.sum()) offspring = np.zeros(n_e, int) for j, k in enumerate(u.unit[Es]): offspring[j] = int(par_e.get(k, 0)) chain_units = ["Southern Cooks", "Society Islands", "Tuamotus", "Mangareva (Gambier)"] res["family_tree"] = dict( documented_links=int(par.sum()), east_units=n_e, east_units_with_documented_daughter=int((offspring > 0).sum()), documented_coloniser_share_lower_bound=float((offspring > 0).mean()), documented_coloniser_share_jeffreys95=jeffreys(int((offspring > 0).sum()), n_e), daughters_per_documented_parent={k: int(v) for k, v in par.items()}, units_with_unknown_parent=int(n_e - par.sum()), offspring_mean=float(offspring.mean()), offspring_variance=float(offspring.var(ddof=1)), early_chain_units_documented_daughters={k: int(par_e.get(k, 0)) for k in chain_units}, early_chain_units_mean_daughters=float(np.mean([par_e.get(k, 0) for k in chain_units])), note="lower bound: the source of %d East units is not established, so their parents are uncounted. Over a finished wave " "on a finite set of targets the mean number of daughters per unit is below one by construction; the direct count " "that matters is the number of daughters of the EARLY units" % int(n_e - par.sum())) # ------------------------------------------------------------------ abandonment ab = (u.status == "abandoned").values k, m = int((E & ab).sum()), int(Es.sum()) variants = {} variants["conventional units (main)"] = (k, m) canon = E & ab & (u.canonical_mystery == 1).values # strict: only canonical mystery islands with firm evidence (Howland recoded to traces), non-canonical dropped from the sample strict_ab = canon & (u.unit != "Howland").values variants["canonical list only, Howland dropped"] = (int(strict_ab.sum()), int((E & pers).sum() + strict_ab.sum())) variants["traces counted as abandoned settlements"] = (k + int((E & (u.status == "traces")).sum()), m + int((E & (u.status == "traces")).sum())) variants["all settled units incl. West Polynesia"] = (int(ab.sum()), int(settled.sum())) for L in (150, 300, 500): lab = single_linkage(D, np.where(Es)[0], L) cl = pd.DataFrame(dict(i=list(lab.keys()), c=list(lab.values()))) cl["pers"] = pers[cl.i.values] g = cl.groupby("c").pers.any() variants[f"clusters, single linkage {L} km (East)"] = (int((~g).sum()), int(len(g))) variants["doubtful settlements (Howland, Necker, Auckland Islands) dropped"] = (k - int((E & ab & doubtful).sum()), m - int((E & ab & doubtful).sum())) hab_only = (u.evidence == "habitation").values variants["only abandoned units with documented habitation sites"] = (int((E & ab & hab_only).sum()), int((E & pers).sum() + (E & ab & hab_only).sum())) variants["undated West units moved into the pulse domain (core sources only)"] = (int((dom_core & ab).sum()), int(dom_core.sum())) res["abandonment"] = dict(fractions={kk: dict(abandoned=v[0], settled=v[1], fraction=v[0] / v[1], jeffreys95=jeffreys(v[0], v[1])) for kk, v in variants.items()}) area = u.land_km2.values res["abandonment"]["area_fraction_abandoned_east"] = float(area[E & ab].sum() / area[Es].sum()) res["abandonment"]["area_fraction_abandoned_east_excl_NZ"] = float( area[E & ab].sum() / area[Es & (u.unit != "New Zealand").values].sum()) # the Auckland Islands 'settlement' is one site on Enderby Island (7.69 km2, Wikidata), not the 570 km2 group area_site = area.copy() area_site[ix["Auckland Islands"]] = 7.69 res["abandonment"]["area_fraction_abandoned_east_enderby_only"] = float(area_site[E & ab].sum() / area_site[Es].sum()) res["abandonment"]["area_fraction_abandoned_east_excl_NZ_enderby_only"] = float( area_site[E & ab].sum() / area_site[Es & (u.unit != "New Zealand").values].sum()) res["abandonment"]["new_zealand_share_of_settled_land"] = float(area[ix["New Zealand"]] / area[Es].sum()) res["abandonment"]["median_area_km2"] = dict(abandoned=float(np.median(area[E & ab])), persisted=float(np.median(area[E & pers]))) res["abandonment"]["median_isolation_big_km"] = dict(abandoned=float(np.median(iso_big[E & ab])), persisted=float(np.median(iso_big[E & pers]))) from scipy.stats import mannwhitneyu res["abandonment"]["mannwhitney_p_area"] = float(mannwhitneyu(area[E & ab], area[E & pers]).pvalue) res["abandonment"]["mannwhitney_p_isolation"] = float(mannwhitneyu(iso_big[E & ab], iso_big[E & pers]).pvalue) def fit(mask, iso, label, ridge=0.0, area_=area): X = np.column_stack([np.log10(area_[mask]), np.log10(iso[mask])]) y = ab[mask].astype(float) b, cov, ll = logistic_fit(X, y, ridge=ridge) b0, _, ll0 = logistic_fit(np.zeros((len(y), 0)), y) ba, _, lla = logistic_fit(X[:, :1], y) # bootstrap bs = [] r2 = np.random.default_rng(SEED + 1) for _ in range(2000): j = r2.integers(0, len(y), len(y)) if y[j].min() == y[j].max(): continue bb, _, _ = logistic_fit(X[j], y[j], ridge=max(ridge, 0.05)) bs.append(bb) bs = np.array(bs) # leave-one-out pl = np.zeros(len(y)) for j in range(len(y)): keep = np.arange(len(y)) != j bb, _, _ = logistic_fit(X[keep], y[keep], ridge=max(ridge, 0.05)) pl[j] = 1 / (1 + np.exp(-(bb[0] + X[j] @ bb[1:]))) pos, neg = pl[y == 1], pl[y == 0] auc = float(((pos[:, None] > neg[None, :]).sum() + 0.5 * (pos[:, None] == neg[None, :]).sum()) / (len(pos) * len(neg))) def pred(a_km2, iso_km, coef=b): z = coef[0] + coef[1] * np.log10(a_km2) + coef[2] * np.log10(iso_km) return float(1 / (1 + np.exp(-z))) def pred_ci(a_km2, iso_km): z = bs[:, 0] + bs[:, 1] * np.log10(a_km2) + bs[:, 2] * np.log10(iso_km) pp = 1 / (1 + np.exp(-z)) return [float(np.percentile(pp, 2.5)), float(np.percentile(pp, 97.5))] cases = {"10 km2 at 300 km": (10, 300), "10 km2 at 1500 km": (10, 1500), "100 km2 at 300 km": (100, 300), "100 km2 at 1500 km": (100, 1500), "1000 km2 at 1500 km": (1000, 1500), "1000 km2 at 300 km": (1000, 300)} return dict(label=label, n=int(len(y)), n_abandoned=int(y.sum()), coef=dict(intercept=float(b[0]), log10_area=float(b[1]), log10_isolation=float(b[2])), se=dict(intercept=float(np.sqrt(cov[0, 0])), log10_area=float(np.sqrt(cov[1, 1])), log10_isolation=float(np.sqrt(cov[2, 2]))), boot95=dict(log10_area=[float(np.percentile(bs[:, 1], 2.5)), float(np.percentile(bs[:, 1], 97.5))], log10_isolation=[float(np.percentile(bs[:, 2], 2.5)), float(np.percentile(bs[:, 2], 97.5))]), loglik=float(ll), loglik_null=float(ll0), loglik_area_only=float(lla), lr_chi2_vs_null=float(2 * (ll - ll0)), lr_chi2_isolation_given_area=float(2 * (ll - lla)), loo_auc=auc, loo_brier=float(np.mean((pl - y) ** 2)), null_brier=float(np.mean((y.mean() - y) ** 2)), predictions={kk: dict(p=pred(*v), boot95=pred_ci(*v)) for kk, v in cases.items()}, area_at_p50_for_isolation={str(d): float(10 ** (-(b[0] + b[2] * np.log10(d)) / b[1])) for d in (300, 1000, 2000)}, ), b fits = {} fits["east_main"], b_main = fit(Es, iso_big, "East Polynesia settled units; isolation = distance to nearest other unit of at least 100 km2") fits["east_iso_persisted"], _ = fit(Es, iso_pers, "East; isolation = distance to nearest persisted unit") fits["all_settled"], _ = fit(settled, iso_big, "West + East settled units") m_strict = (E & pers) | strict_ab fits["east_canonical_only"], _ = fit(m_strict, iso_big, "East; abandoned = canonical list without Howland") m_notsub = Es & (u.type != "subantarctic").values fits["east_no_subantarctic"], _ = fit(m_notsub, iso_big, "East; Auckland Islands removed") fits["east_no_doubtful"], _ = fit(Es & ~doubtful, iso_big, "East; Howland, Necker and Auckland Islands removed") res["abandonment"]["logistic"] = fits res["abandonment"]["_coef_main"] = [float(x) for x in b_main] res["abandonment"]["unit_level_caveat"] = ( "the unit table merges persisted small islands into archipelago units while every abandoned island is its own unit, so " "'area' at unit level is the land of a whole island cluster and absorbs the effect of having neighbours; the fitted " "chances of abandonment therefore apply to stand-alone islands and to whole archipelagos, not to an island inside one") # island-level analysis (revision) res["abandonment"]["island_level"], isl = island_level(u, a, D) # settlement lifetimes of abandoned units that have both dates t = sample_dates(u, "S", rng, 20000) life = {} for i, r in u.iterrows(): if r.status == "abandoned" and np.isfinite(r.last_occ) and np.isfinite(r.S_lo): L = r.last_occ - t[:, i] life[r.unit] = [float(np.percentile(L, p)) for p in (16, 50, 84)] res["abandonment"]["lifetime_years_p16_p50_p84"] = life res["abandonment"]["lifetime_median_of_medians"] = float(np.median([v[1] for v in life.values()])) usable = [k_ for k_ in life if (u.S_lo[ix[k_]], u.S_hi[ix[k_]]) not in GENERIC_WINDOWS] res["abandonment"]["lifetime_usable_dates"] = dict( units=usable, n=len(usable), median_range_yr=[float(min(life[k_][1] for k_ in usable)), float(max(life[k_][1] for k_ in usable))], p16_p84_envelope_yr=[float(min(life[k_][0] for k_ in usable)), float(max(life[k_][2] for k_ in usable))], dropped_generic_window=[k_ for k_ in life if k_ not in usable], note="lifetime = last dated occupation minus a first-settlement window; indicative only") # crude failure hazard: failures per settlement-year of exposure, by island size. # Exposure runs from the median first-settlement date (chronology S) to the last dated occupation, or to 1770 CE # (about the time of European contact) for units still inhabited. Units without dates get the median exposure of # their own outcome class; that assumption is recorded in the output. tmed = np.nanmedian(t, axis=0) expo = np.where(ab, u.last_occ.values - tmed, 1770.0 - tmed) haz = {} for label, msk in (("all East units", Es), ("under 100 km2", Es & (area < 100)), ("100 km2 or more", Es & (area >= 100))): e = expo.copy() for cls in (ab, pers): sel = msk & cls fill = np.nanmedian(expo[Es & cls]) e[sel & ~np.isfinite(e)] = fill kf, ex = int((msk & ab).sum()), float(e[msk].sum()) from scipy.stats import gamma lo_, hi_ = gamma.ppf([0.025, 0.975], kf + 0.5) / ex haz[label] = dict(units=int(msk.sum()), failures=kf, exposure_settlement_years=ex, hazard_per_yr=kf / ex, hazard_per_yr_95=[float(lo_), float(hi_)], mean_time_to_failure_yr=ex / max(kf, 1), units_with_imputed_exposure=int((msk & ~np.isfinite(expo)).sum())) res["abandonment"]["hazard"] = haz res["abandonment"]["relaunch_to_lifetime_ratio"] = dict( relaunch_delay_yr=res["timing"]["genetic_chain"]["relaunch_delay_median_yr"], failed_settlement_lifetime_yr=res["abandonment"]["lifetime_median_of_medians"], ratio=res["timing"]["genetic_chain"]["relaunch_delay_median_yr"] / res["abandonment"]["lifetime_median_of_medians"], note="indicative only: eight poorly dated lifetimes divided by a delay taken from four point dates of the other chronology") # ------------------------------------------------------------------ reach: share of targets settled, by definition X = (u.region == "X").values later_s = (u.later_strict == 1).values later_l = (u.later_lenient == 1).values harsh = u.type.isin(["subantarctic", "rock"]).values in_reach = E & (nearest_settled <= D_net) def share(den, label, outcome_dependent=False, note=""): kk, nn = int((den & settled).sum()), int(den.sum()) k2 = int((den & settled & ~doubtful).sum()) return dict(definition=label, targets=nn, settled=kk, fraction=kk / nn, jeffreys95=jeffreys(kk, nn), inhabited_at_contact=int((den & pers).sum()), settled_without_doubtful=k2, unsettled=sorted(u.unit[den & ~settled]), denominator_depends_on_outcome=outcome_dependent, note=note) defs = [ share(in_reach, "every land unit of East Polynesia in reach"), share(in_reach & ~harsh, "excluding rocks and subantarctic islands"), share(in_reach & ~harsh & (area >= 10), "land of 10 km2 or more, excluding rocks and subantarctic islands"), share(in_reach & later_l, "carried a civilian settlement for about 20 years or more after European contact (lenient later-habitability)"), share(in_reach & later_s, "carries a self-supporting community founded or continued after European contact (strict later-habitability)"), share(in_reach & (settled | later_s), "first-pass definition: settled prehistorically OR strict later-habitability", True, "29 of the 30 'trials' are in the denominator because they succeeded; no interval is meaningful"), ] res["reach"] = dict( definition="'in reach' = within the network range of a settled unit (all 50 East units qualify); see reach_from_sources for a " "version that does not use the settled set at all", by_definition=defs, beyond_network=dict(later_habitable_units={k_: float(nearest_settled[ix[k_]]) for k_ in u.unit[X & later_s]}, settled=int((X & later_s & settled).sum()), n=int((X & later_s).sum())), ) # outcome-independent reach: units reachable from West Polynesia by hops of at most R over ANY land (settled or not) anyland = np.ones(n, bool) bn_any, _ = bottleneck_from_sources(D, anyland, sources) rfs = {} for R in (2100, 2600, 3500, 4000): inr = (E | X) & (bn_any <= R) rfs[str(R)] = dict(land_units=int(inr.sum()), settled=int((inr & settled).sum()), strict_later_habitable=int((inr & later_s).sum()), strict_later_habitable_settled=int((inr & later_s & settled).sum()), strict_later_habitable_unsettled=sorted(u.unit[inr & later_s & ~settled]), lenient_later_habitable=int((inr & later_l).sum()), lenient_later_habitable_settled=int((inr & later_l & settled).sum())) res["reach"]["reach_from_sources"] = rfs # what predicts WHETHER a land unit was settled? (all 50 East units; predictors do not depend on the outcome) Xs = np.column_stack([np.log10(area[E]), harsh[E].astype(float)]) ys = settled[E].astype(float) bset, covset, llset = logistic_fit(Xs, ys) _, _, ll0set = logistic_fit(np.zeros((len(ys), 0)), ys) pl = np.zeros(len(ys)) for j in range(len(ys)): keep = np.arange(len(ys)) != j bb, _, _ = logistic_fit(Xs[keep], ys[keep], ridge=0.05) pl[j] = _sigmoid(bb[0] + Xs[j] @ bb[1:]) pos, neg = pl[ys == 1], pl[ys == 0] cls = [(0, 1), (1, 10), (10, 100), (100, 1e9)] res["reach"]["settlement_vs_size"] = dict( n=int(len(ys)), settled=int(ys.sum()), coef=dict(intercept=float(bset[0]), log10_area=float(bset[1]), subantarctic_or_rock=float(bset[2])), se=dict(log10_area=float(np.sqrt(covset[1, 1])), subantarctic_or_rock=float(np.sqrt(covset[2, 2]))), lr_chi2_vs_null=float(2 * (llset - ll0set)), loo_auc=float(((pos[:, None] > neg[None, :]).sum() + 0.5 * (pos[:, None] == neg[None, :]).sum()) / (len(pos) * len(neg))), by_area_class=[dict(area_km2=f"{lo}-{hi if hi < 1e8 else 'inf'}", units=int((E & (area >= lo) & (area < hi)).sum()), settled=int((Es & (area >= lo) & (area < hi)).sum()), inhabited_at_contact=int((E & pers & (area >= lo) & (area < hi)).sum())) for lo, hi in cls], unsettled_units_km2={k_: float(area[ix[k_]]) for k_ in u.unit[E & ~settled]}, median_area_unsettled_km2=float(np.median(area[E & ~settled])), median_area_settled_km2=float(np.median(area[Es]))) # kept for backward compatibility with the first-pass key numbers (flagged as outcome-dependent in the report) kf, nf = defs[-1]["settled"], defs[-1]["targets"] res["reach"]["fill_fraction"] = dict(k=kf, n=nf, fraction=kf / nf, outcome_dependent=True) res["reach"]["settled_fraction_of_all_islands_in_range"] = dict( k=defs[0]["settled"], n=defs[0]["targets"], fraction=defs[0]["fraction"], persisted_fraction=defs[0]["inhabited_at_contact"] / defs[0]["targets"]) # ------------------------------------------------------------------ percolation fit # Node set A (first pass): the settled units plus Lord Howe. It is CONDITIONED ON THE OUTCOME. # Node set B (revision): the same simulation, but every unsettled East land unit enters the likelihood through a # habitability term that depends only on land area and on being subantarctic or a rock. nodes = np.where(settled | (u.unit == "Lord Howe").values)[0] Dn = D[np.ix_(nodes, nodes)] nm = [names[i] for i in nodes] jx = {k: j for j, k in enumerate(nm)} settled0 = np.array([W[i] for i in nodes]) must = np.array([j for j, i in enumerate(nodes) if Es[i]]) mustnot = np.array([jx["Lord Howe"]]) others = np.where(E & ~settled & (u.unit != "Lord Howe").values)[0] extra = dict(D_extra=D[np.ix_(nodes, others)], x_set=np.log10(area[nodes[must]]), z_set=harsh[nodes[must]].astype(float), x_un=np.log10(area[np.concatenate([nodes[mustnot], others])]), z_un=harsh[np.concatenate([nodes[mustnot], others])].astype(float)) p_grid = np.round(np.arange(0.05, 1.0001, 0.05), 2) s_grid = np.array([0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) perc = {} rows = [] i35 = int(np.where(np.isclose(p_grid, 0.35))[0][0]) for R in (2100, 2600, 3500): cover, exact, frac, likB, parB = percolation_mc(Dn, settled0, jx["Samoa"], must, mustnot, p_grid, s_grid, R, rng, extra=extra) deg = ((Dn <= R).sum(1) - 1) allE = np.where(E | W)[0] deg_all = ((D[np.ix_(allE, allE)] <= R).sum(1) - 1) zc = 4.512 # critical mean degree of a 2-D random geometric graph (continuum percolation of discs) out = dict(range_km=R, mean_degree=float(deg.mean()), min_degree=int(deg.min()), p_c_site_rgg2d=float(min(1.0, zc / deg.mean())), mean_degree_all_land_units=float(deg_all.mean()), p_c_all_land_units=float(min(1.0, zc / deg_all.mean())), node_set_note="degrees on the settled-plus-Lord-Howe graph (selected on the outcome) and on all land units") for label, Lk in (("exact_pattern", exact), ("coverage_only", cover)): post = Lk / Lk.sum() if Lk.sum() > 0 else Lk mp = post.sum(1) ms = post.sum(0) cdf = np.cumsum(mp) out[label] = dict( max_likelihood=float(Lk.max()), p_mle=float(p_grid[np.unravel_index(Lk.argmax(), Lk.shape)[0]]), s_mle=float(s_grid[np.unravel_index(Lk.argmax(), Lk.shape)[1]]), p_post_mean=float((mp * p_grid).sum()) if Lk.sum() > 0 else None, p_post_q05=float(p_grid[np.searchsorted(cdf, 0.05)]) if Lk.sum() > 0 else None, p_post_q50=float(p_grid[np.searchsorted(cdf, 0.50)]) if Lk.sum() > 0 else None, s_post_mean=float((ms * s_grid).sum()) if Lk.sum() > 0 else None, post_prob_p_below_0p35=float(mp[p_grid <= 0.35].sum()) if Lk.sum() > 0 else None, lik_ratio_p0p35_vs_best=float(Lk[np.isclose(p_grid, 0.35)].max() / Lk.max()) if Lk.max() > 0 else None, strict_landis_s1_max_lik=float(Lk[:, -1].max()), ) # strict Landis rule (s = 1), coverage only if label == "coverage_only": l1 = Lk[:, -1] c1 = np.cumsum(l1 / l1.sum()) out["strict_landis_s1_coverage"] = dict(p_post_q05=float(p_grid[np.searchsorted(c1, 0.05)]), p_post_mean=float((l1 / l1.sum() * p_grid).sum()), lik_at_p0p35=float(l1[np.isclose(p_grid, 0.35)][0]), lik_at_p1=float(l1[-1])) # direct, non-asymptotic statement: how far does a wave with p = 0.35 get on this graph? out["wave_at_p0p35"] = dict( strict_landis_s1=dict(mean_fraction_of_units_reached=float(frac[i35, -1]), prob_all_reached=float(cover[i35, -1])), s0p5=dict(mean_fraction_of_units_reached=float(frac[i35, list(s_grid).index(0.5)]), prob_all_reached=float(cover[i35, list(s_grid).index(0.5)]))) out["smallest_p_reaching_90pct_of_units_on_average_s1"] = float(p_grid[np.argmax(frac[:, -1] >= 0.9)]) if (frac[:, -1] >= 0.9).any() else None # model B: profile likelihood of p (maximised over s and the habitability parameters) profB = likB.max(1) jb = int(profB.argmax()) sb = int(likB[jb].argmax()) sup = p_grid[profB >= profB.max() * np.exp(-2.0)] th = parB[jb, sb] out["all_land_units_with_habitability"] = dict( n_settled=int(len(must)), n_unsettled=int(1 + len(others)), max_likelihood=float(profB.max()), p_at_max=float(p_grid[jb]), s_at_max=float(s_grid[sb]), habitability_at_max=dict(intercept=float(th[0]), per_dex_area=float(th[1]), subantarctic_or_rock=float(th[2]), h_at_1km2=float(_sigmoid(th[0])), h_at_10km2=float(_sigmoid(th[0] + th[1])), h_at_100km2=float(_sigmoid(th[0] + 2 * th[1]))), profile_lik_ratio_p0p35_vs_best=float(profB[i35] / profB.max()), p_support_within_e2=[float(sup.min()), float(sup.max())], profile={str(p_): float(v_ / profB.max()) for p_, v_ in zip(p_grid, profB)}, strict_landis_s1_max_lik_ratio=float(likB[:, -1].max() / profB.max()), flat_prior_mean_p=float((profB / profB.sum() * p_grid).sum()), ) perc[str(R)] = out for a_, p in enumerate(p_grid): for b_, s in enumerate(s_grid): rows.append(dict(range_km=R, p=p, s=s, lik_coverage=cover[a_, b_], lik_exact=exact[a_, b_], mean_fraction_reached=frac[a_, b_], lik_all_land_habitability=likB[a_, b_])) pd.DataFrame(rows).to_csv(os.path.join(DATA, "pacific_percolation_likelihood_grid.csv"), index=False) # combined posterior over the three ranges (equal weight), exact pattern g = pd.DataFrame(rows) comb = g.groupby("p").lik_exact.sum() comb = comb / comb.sum() cdf = comb.cumsum() perc["combined_exact"] = dict(p_post_mean=float((comb * comb.index).sum()), p_post_q05=float(cdf.index[np.searchsorted(cdf.values, 0.05)]), p_post_q50=float(cdf.index[np.searchsorted(cdf.values, 0.5)]), post_prob_p_below_0p35=float(comb[comb.index <= 0.35].sum()), flat_prior_mean_on_this_grid=float(p_grid.mean()), flat_prior_prob_p_below_0p35=float((p_grid <= 0.35).mean()), warning="prior-dominated and conditioned on a node set made of the settled units plus Lord Howe") # gateways: settled East units with very few settled neighbours in range gate = {} for R in (2100, 2600): gate[str(R)] = {r.unit: [names[j] for j in np.where(settled & (Dm[i] <= R))[0]] for i, r in u.iterrows() if Es[i] and (settled & (Dm[i] <= R)).sum() <= 3} perc["few_parent_units"] = gate res["percolation"] = perc u.to_csv(os.path.join(DATA, "pacific_units_derived.csv"), index=False) a.to_csv(os.path.join(DATA, "pacific_anchor_islands_with_coordinates.csv"), index=False) dump_json(res, "pacific_results.json") if verbose: print(json.dumps({k: res[k] for k in ("counts", "hops")}, indent=1, default=float)[:6000]) return res, u, D, isl if __name__ == "__main__": run()