#!/usr/bin/env python3 """Single entry point for the expansion_analogs study. python3 run_all.py # regenerate every table, number and figure python3 run_all.py --download # force a fresh download of the four public raw files first python3 run_all.py --no-model # skip the re-run of the project's integrated model (saves about 25 s) python3 run_all.py --refresh-geo # re-query Wikidata for island coordinates (slow: the endpoint is rate-limited) python3 run_all.py --refresh-islands # re-download the Wikipedia infoboxes behind the island-level table (about 3 minutes) Inputs code/curated_tables.py hand-entered tables with a source key and verification grade on every row data/wikidata_geo_cache.json island coordinates from Wikidata (cached; --refresh-geo rebuilds it) data/island_infobox_cache.json coordinates and land areas parsed from Wikipedia infoboxes (cached; --refresh-islands rebuilds it) raw downloads fermi_raw/expansion_analogs/ (fetched automatically if missing) Outputs data/curated/*.csv, data/*.json, data/*.csv, figures/*.png, data/key_numbers.json """ import argparse import json import os import sys import time import warnings HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, HERE) sys.dont_write_bytecode = True import common # noqa: E402 def refresh_geo(): """Re-query Wikidata for the coordinates (property P625) of every anchor island, using the item ids stored in data/wikidata_geo_cache.json (the ids were found once from the English Wikipedia titles used as geo_key). One SPARQL request; the public endpoint is rate-limited, so failures are retried after a pause.""" import urllib.parse import urllib.request import curated_tables as ct keys = sorted({r[2] for r in ct.PACIFIC_ANCHORS if r[2]}) path = os.path.join(common.DATA, "wikidata_geo_cache.json") old = json.load(open(path)) qid = {k: old[k]["qid"] for k in keys} q = "SELECT ?item ?c WHERE { VALUES ?item { %s } ?item wdt:P625 ?c . }" % " ".join("wd:" + v for v in sorted(set(qid.values()))) url = "https://query.wikidata.org/sparql?format=json&query=" + urllib.parse.quote(q) req = urllib.request.Request(url, headers={"User-Agent": "expansion_analogs research script (python-urllib)"}) for attempt in range(5): try: with urllib.request.urlopen(req, timeout=120) as r: rows = json.load(r)["results"]["bindings"] break except Exception as e: # rate limit: wait and retry print(" wikidata retry:", e) time.sleep(65) else: raise RuntimeError("Wikidata query failed; the cached coordinates were left untouched") got = {} for b in rows: lon, lat = b["c"]["value"].replace("Point(", "").replace(")", "").split() got.setdefault(b["item"]["value"].rsplit("/", 1)[-1], []).append([float(lat), float(lon)]) missing = [k for k in keys if qid[k] not in got] if missing: raise RuntimeError("no coordinates returned for: %s" % missing) moved = {} for k in keys: cands = got[qid[k]] o = old[k]["wd_coord"] # an item can carry several coordinate statements: keep the one closest to the cached value new = min(cands, key=lambda c: abs(c[0] - o[0]) + abs(((c[1] - o[1] + 180) % 360) - 180)) if max(abs(new[0] - o[0]), abs(((new[1] - o[1] + 180) % 360) - 180)) > 0.05: moved[k] = (o, new) old[k]["wd_coord"] = new json.dump(old, open(path, "w"), indent=1, ensure_ascii=False) print(" coordinates refreshed for %d islands; %d moved by more than 0.05 degrees %s" % (len(keys), len(moved), moved)) def key_numbers(pac, norse, inv, syn): """Flat table of the numbers quoted in the report.""" t, ab, pe, hp = pac["timing"], pac["abandonment"], pac["percolation"], pac["hops"] f = ab["fractions"] k = { "pacific_units": pac["counts"]["n_units"], "pacific_anchor_islands": pac["counts"]["n_anchor_islands"], "east_settled": pac["counts"]["east_settled"], "east_persisted": pac["counts"]["east_persisted"], "east_abandoned": pac["counts"]["east_abandoned"], "verification_grades": pac["counts"]["verification_grades"], "abandoned_evidence": pac["counts"]["abandoned_evidence"], "doubtful_settlements": pac["counts"]["doubtful_settlements"], "hop_nearest_median_km": hp["nearest_settled_neighbour_km"]["median"], "hop_nearest_max_km": hp["nearest_settled_neighbour_km"]["max"], "hop_nearest_quartiles_km": [hp["nearest_settled_neighbour_km"]["q25"], hp["nearest_settled_neighbour_km"]["q75"]], "hop_evidence_n": hp["evidence_based_hops_km"]["n"], "hop_evidence_median_km": hp["evidence_based_hops_km"]["median"], "hop_evidence_quartiles_km": [hp["evidence_based_hops_km"]["q25"], hp["evidence_based_hops_km"]["q75"]], "hop_evidence_max_km": hp["evidence_based_hops_km"]["max"], "network_bottleneck_km": hp["network_bottleneck_km"], "corners": hp["corners"], "rapa_nui_to_south_america_km": hp["rapa_nui_to_south_america_km"], "attested_range_summary": hp["attested_range_summary"], "core_sources_sensitivity": hp["core_sources_sensitivity"], "frac_nearest_hops_below_half_range": hp["frac_nearest_hops_below_half_bottleneck"], "frac_evidence_hops_below_half_range": hp["frac_evidence_hops_below_half_bottleneck"], } for c in ("S", "A"): for q in ("first_east_settlement", "pause_years", "pulse_10_90_years", "doubling_time_first_to_half_yr", "R_per_30yr_generation_first_to_half", "doubling_time_burst_yr", "R_per_30yr_generation_burst", "new_units_after_1350", "units_by_1350"): k[f"{c}_{q}_p2.5_p16_p50_p84_p97.5"] = t[c][q] k[f"{c}_relaunch_delay_east_parents_pooled"] = t[c]["relaunch_delay_east_parents_pooled"] k[f"{c}_relaunch_delay_frac_under_100yr"] = t[c]["relaunch_delay_east_parents_frac_under_100yr"] k[f"{c}_relaunch_delay_by_parent_median"] = t[c]["relaunch_delay_by_parent_median"] k["growth_factor_note"] = t["growth_note"] k["growth_sensitivity"] = t["growth_sensitivity"] k["generic_window_units"] = t["generic_window_units"] k["genetic_chain"] = t["genetic_chain"] k["family_tree"] = pac["family_tree"] for c in ("S", "A"): k[f"{c}_growth_factor_per_relaunch_delay"] = t[c]["R_per_relaunch_delay_first_to_half"] k[f"{c}_late_new_units_per_existing_unit"] = t[c]["late_new_units_per_existing_unit"] k["abandonment_fractions"] = f k["abandonment_area_fraction"] = ab["area_fraction_abandoned_east"] k["abandonment_area_fraction_excl_NZ"] = ab["area_fraction_abandoned_east_excl_NZ"] k["abandonment_area_fraction_enderby_only"] = ab["area_fraction_abandoned_east_enderby_only"] k["abandonment_area_fraction_excl_NZ_enderby_only"] = ab["area_fraction_abandoned_east_excl_NZ_enderby_only"] k["new_zealand_share_of_settled_land"] = ab["new_zealand_share_of_settled_land"] k["abandonment_median_area_km2"] = ab["median_area_km2"] k["abandonment_mannwhitney_p_area"] = ab["mannwhitney_p_area"] k["abandonment_mannwhitney_p_isolation"] = ab["mannwhitney_p_isolation"] k["abandonment_logistic_unit_level"] = {kk: dict(coef=v["coef"], se=v["se"], boot95=v["boot95"], lr_chi2_vs_null=v["lr_chi2_vs_null"], lr_chi2_isolation_given_area=v["lr_chi2_isolation_given_area"], loo_auc=v["loo_auc"], n=v["n"], n_abandoned=v["n_abandoned"]) for kk, v in ab["logistic"].items()} k["abandonment_predictions_unit_level_main"] = ab["logistic"]["east_main"]["predictions"] k["abandonment_unit_level_caveat"] = ab["unit_level_caveat"] k["abandonment_island_level"] = ab["island_level"] k["abandonment_lifetimes"] = ab["lifetime_years_p16_p50_p84"] k["abandonment_lifetime_median_of_medians"] = ab["lifetime_median_of_medians"] k["abandonment_lifetime_usable_dates"] = ab["lifetime_usable_dates"] k["abandonment_hazard"] = ab["hazard"] k["relaunch_to_lifetime_ratio"] = ab["relaunch_to_lifetime_ratio"] k["reach"] = pac["reach"] k["percolation"] = {R: dict(mean_degree=v["mean_degree"], p_c=v["p_c_site_rgg2d"], mean_degree_all_land_units=v["mean_degree_all_land_units"], p_c_all_land_units=v["p_c_all_land_units"], exact=v["exact_pattern"], coverage=v["coverage_only"], strict_landis=v.get("strict_landis_s1_coverage"), wave_at_p0p35=v["wave_at_p0p35"], all_land_units_with_habitability=v["all_land_units_with_habitability"]) for R, v in pe.items() if R.isdigit()} k["percolation_combined_first_pass_node_set"] = pe["combined_exact"] k["norse"] = {x: norse[x] for x in ("n_frontier_units", "abandoned", "jeffreys95", "empty_land_units", "empty_land_abandoned", "hop_km", "relaunch_lag_years_p16_p50_p84", "lifetimes_years", "wave")} k["invasions"] = dict(propagule={x: inv["propagule"][x] for x in ("n_events", "n_species", "overall_fraction", "overall_jeffreys95", "bins", "logistic_all", "predicted_establishment", "species_concentration", "literature_note")}, gavia=inv["gavia"], cassey=inv["cassey"], p3k14c=inv["p3k14c"]) k["synthesis"] = {x: syn[x] for x in syn if x != "model_rerun"} if "model_rerun" in syn and "error" not in syn["model_rerun"]: k["model_rerun_sensitivity_test"] = { name: {s: dict(posterior=v["schools"][s]["posterior"], bayes_factor_base_over_this=v["schools"][s].get("bayes_factor_base_over_this")) for s in v["schools"]} for name, v in syn["model_rerun"].items() if name != "class_labels"} return k def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--download", action="store_true", help="force re-download of the raw public data files") ap.add_argument("--no-model", action="store_true", help="skip the integrated-model re-run") ap.add_argument("--refresh-geo", action="store_true", help="re-query Wikidata for island coordinates") ap.add_argument("--refresh-islands", action="store_true", help="re-download and re-parse the Wikipedia infoboxes of the island table") ap.add_argument("--quiet", action="store_true") a = ap.parse_args() v = not a.quiet t0 = time.time() for d in (common.DATA, common.CUR, common.FIG): os.makedirs(d, exist_ok=True) import curated_tables curated_tables.write_all(common.CUR) print("[1/7] curated tables written to", common.CUR) for key in common.DOWNLOADS: try: p = common.fetch(key, force=a.download) print(" raw file %-12s %s (%.1f MB)" % (key, p, os.path.getsize(p) / 1e6)) except Exception as e: print(" raw file %-12s DOWNLOAD FAILED: %s" % (key, e)) if key != "p3k14c": raise if a.refresh_geo: refresh_geo() import island_sources if a.refresh_islands or not os.path.exists(island_sources.CACHE): island_sources.build_cache(force=a.refresh_islands) print(" island infobox cache rebuilt (%s)" % island_sources.CACHE) print("[2/7] raw data present (%.0f s)" % (time.time() - t0)) import pacific, norse, invasions, synthesis, figures with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) pac, u, D, isl = pacific.run(verbose=False) print("[3/7] Pacific analysis done (%.0f s)" % (time.time() - t0)) nor = norse.run(verbose=v) print("[4/7] Norse analysis done (%.0f s)" % (time.time() - t0)) inv, tab, coefs = invasions.run(verbose=v) print("[5/7] invasion base rates done (%.0f s)" % (time.time() - t0)) syn = synthesis.run(pac, nor, inv, verbose=v, with_model=not a.no_model) print("[6/7] synthesis done (%.0f s)" % (time.time() - t0)) with warnings.catch_warnings(): warnings.simplefilter("ignore") figures.run(u, pac, nor, inv, tab, coefs, syn, isl) kn = key_numbers(pac, nor, inv, syn) common.dump_json(kn, "key_numbers.json") print("[7/7] figures written to %s; key numbers in data/key_numbers.json (%.0f s total)" % (common.FIG, time.time() - t0)) if __name__ == "__main__": main()