"""Single entry point of the planet_lifetime study. Regenerates every number and figure. python run_all.py # downloads missing raw tables (VizieR TAP), then runs everything (see report for run time) python run_all.py --refresh # force re-download of all raw tables python run_all.py --quick # small Monte Carlo sizes, for testing only (numbers will differ slightly) Outputs: data/*.csv, data/results_summary.json, figures/fig*.png """ import json import os import sys import time import numpy as np import analysis as an import fetch_data import figures import population as pop import stellar import thermal import verify from common import DATA_DIR, AGE_SUN def _clean(o): if isinstance(o, dict): return {str(k): _clean(v) for k, v in o.items()} if isinstance(o, (list, tuple)): return [_clean(v) for v in o] if isinstance(o, (np.bool_,)): return bool(o) if isinstance(o, (np.floating, np.integer)): return float(o) if isinstance(o, np.ndarray): return o.tolist() return o def main(): quick = "--quick" in sys.argv t0 = time.time() def log(msg): print(f"[{time.time() - t0:6.0f}s] {msg}", flush=True) log("stage 0: raw data") fetch_data.fetch_all(refresh="--refresh" in sys.argv, verbose=not quick) log("stage 1: abundances") A = an.abundance_stage(n_boot=60 if quick else 300) S, ab_table, ab_summary, emp = A["S"], A["table"], A["summary"], A["emp"] print(A["raw_slopes"].round(4).to_string()) print(A["scatter"].round(4).to_string()) print(A["bivariate"].round(4).to_string()) print(A["binned"].round(3).to_string()) print(A["eu_mg"].round(3).to_string()) print("adopted abundance hyper-parameters:", {k: round(v, 4) for k, v in emp.items()}) heat_model, heat_stars = an.heat_budget_stage(emp) print(heat_model[heat_model["scatter_included"].str.startswith("r-process scatter +")].round(3).to_string()) print(heat_stars.round(3).to_string()) log("stage 2: stellar habitable-zone lifetimes") hz, hz_curves = stellar.hz_table() print(hz.round(2).to_string()) log("stage 3: Earth calibration of the thermal models") cal, hist, th_sens = an.earth_stage() print(cal.round(3).to_string()) log("stage 4: baseline population") B = an.baseline_stage(emp, n=1000 if quick else 5000) surv, pct, rk, ref = B["survival"], B["percentiles"], B["ranks"], B["reference"] print(surv.round(4).to_string()) print(rk.round(4).to_string()) print(ref.round(5).to_string()) print(B["peaks"].round(3).to_string()) print(B["tails"].round(3).T.to_string()) log("stage 5: what controls the lifetime") var = an.variance_stage(emp, n=500 if quick else 2000) print(var[var["criterion"] == "degassing>=10%"].round(2).to_string()) log("stage 6: scenarios (tail assumptions, heat budget, thermostat, parameters)") scen, ranges = an.scenario_stage(emp, n=400 if quick else 2000) print(scen[["scenario", "formulation", "criterion", "median_lifetime_solar_epoch", "p95_lifetime_solar_epoch", "P_hab_20Gyr", "rank_n2", "rank_n4", "rank_n6"]].round(4).to_string()) print(ranges.round(4).T.to_string()) log("stage 7: hyper-parameter Monte Carlo") hyper, hsum = an.hyper_stage(emp, n_draws=5 if quick else 60, n=300 if quick else 500) print(hsum.round(4).T.to_string()) shares = an.conditional_shares(hyper) print(shares.round(3).to_string()) log("stage 8: figures") figures.make_all(S, ab_summary, emp) log("stage 9: verification (independent integrator, mixture weights, first-pass regression)") ver = verify.run(emp, n=300 if quick else 600) print(ver.round(5).to_string()) if not ver["pass"].all(): print("WARNING: verification failures:", ver.loc[~ver["pass"], "check"].tolist()) # files written by the first-pass version that no longer exist in this pipeline for stale in ("dissolve_probabilities.csv", "variants.csv"): path = os.path.join(DATA_DIR, stale) if os.path.exists(path): os.remove(path) summary = { "abundance_hyper": emp, "abundance_sample_sizes": ab_summary["n"] | {"th_direct": ab_summary["n_th"], "eu_twins": ab_summary["n_eu_twins"]}, "galah_k_n": {k: {"pass_cuts": v["n_pass_cuts"], "after_clipping": v["n"]} for k, v in ab_summary["fits_k"].items()}, "galah_eu_detections": ab_summary["galah_eu"], "abundance_raw_slopes": A["raw_slopes"].to_dict(orient="records"), "abundance_scatter_robustness": A["scatter"].to_dict(orient="records"), "abundance_bivariate_fits": A["bivariate"].to_dict(orient="records"), "abundance_age_bins_broad": A["binned"].to_dict(orient="records"), "abundance_rock_reference_check": A["rock_refs"].to_dict(orient="records"), "th_eu_birth_slope": {"observed": ab_summary["fit_theu"]["b"], "error": ab_summary["fit_theu"]["e_b"], "one_zone_model": ab_summary["theu_model_slope"]}, "mishenina2022_th_check": A["mishenina"], "eu_mg_observed_ranges": A["eu_mg"].to_dict(orient="records"), "abundance_fit_table": ab_table.to_dict(orient="records"), "initial_heat_budget_model": heat_model.to_dict(orient="records"), "initial_heat_budget_stars": heat_stars.to_dict(orient="records"), "stellar_hz": hz.to_dict(orient="records"), "earth_calibration": cal.to_dict(orient="records"), "survival_baseline_solar_epoch": surv[np.isclose(surv["stellar_age_today_gyr"], AGE_SUN)].to_dict(orient="records"), "survival_baseline_all_epochs": surv.to_dict(orient="records"), "lifetime_percentiles": pct.to_dict(orient="records"), "rank_baseline": rk.to_dict(orient="records"), "rank_reference_fixed_caps": ref.to_dict(orient="records"), "rank_tail_dependence": B["tails"].to_dict(orient="records"), "peak_temperature_baseline": B["peaks"].to_dict(orient="records"), "variance_decomposition": var.to_dict(orient="records"), "scenarios": scen.to_dict(orient="records"), "rank_ranges": ranges.to_dict(orient="records"), "hyper_summary": hsum.to_dict(orient="records"), "conditional_shares": shares.to_dict(orient="records"), "verification": ver.to_dict(orient="records"), "node_ages": pop.NODE_AGES, "target_ages": pop.TARGET_AGES, "quick_mode": quick, "runtime_s": time.time() - t0, } with open(os.path.join(DATA_DIR, "results_summary.json"), "w") as f: json.dump(_clean(summary), f, indent=1) log("done") if __name__ == "__main__": main()