"""Single entry point: regenerates every table, number and figure of the 'convergence' study.

    python3 run_all.py            # fetch missing raw data, rebuild derived tables, run both analyses
    python3 run_all.py --offline  # skip downloads; use the derived CSVs already in data/

Raw downloads go to $FERMI_RAW (default fermi_raw/convergence). Runtime: about 2 minutes
(single thread) once the raw files are present; the downloads themselves (about 30 MB, Europe PMC is
slow) can take several minutes.
"""
import os
import sys
import time

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)


def headline_table(s1, s2):
    """Every number quoted in the first paragraph of the report, taken from the two summaries."""
    T = {t["coding"]: t for t in s1["tests"]}

    def orci(k):
        t = T[k]
        return (f"{t['odds_ratio_per_dex']:.2f} (profile 95% interval {t['odds_ratio_profile_lo95']:.2f}-"
                f"{t['odds_ratio_profile_hi95']:.2f}; permutation p = {t['firth_slope_perm_p']:.3f}; n = {t['n']})")
    il = s1["information_loss"]
    rows = [("Vermeij as published: odds ratio of 'unique' per tenfold age", orci("Vermeij 2006 as published (n=77)")),
            ("  Fisher exact at 543 Ma / 1 Ga", f"{T['Vermeij 2006 as published (n=77)']['fisher_543Ma_table'].replace('/', ' of ').replace(' vs ', ' unique vs ')} "
                                              f"repeated rows predate 543 Ma; p = {T['Vermeij 2006 as published (n=77)']['fisher_543Ma_p']:.5f} / "
                                              f"{T['Vermeij 2006 as published (n=77)']['fisher_1Ga_p']:.4f}"),
            ("  without rows dated >= 3.5 Ga", orci("Vermeij as published, excluding LUCA-era items (>=3.5 Ga)")),
            ("  Phanerozoic only", orci("Vermeij as published, Phanerozoic only (<543 Ma)")),
            ("  Phanerozoic only, without human language", orci("Vermeij as published, Phanerozoic, excluding <5 Ma (language)")),
            ("  Precambrian only", orci("Vermeij as published, Precambrian only (>=543 Ma)"))]
    for lab, c in il["collapsed_coinherited"].items():
        a_, b_ = c['fisher_543Ma_table'].split(" vs ")
        rows.append((f"  {lab}", f"{a_.replace('/', ' of ')} unique vs {b_.replace('/', ' of ')} repeated rows predate 543 Ma, "
                                 f"Fisher p = {c['fisher_543Ma_p']:.3f}; odds ratio "
                                 f"{c['odds_ratio_per_dex']:.2f} ({c['odds_ratio_profile_ci95'][0]:.2f}-{c['odds_ratio_profile_ci95'][1]:.2f})"))
    rows.append(("  step vs gradient log-likelihood (best cut)", f"{il['step_best_cut']['loglik_step']:.1f} vs "
                 f"{il['loglik_gradient_logage']:.1f} (cut between {il['step_best_cut']['between_Ma'][0]:.0f} and "
                 f"{il['step_best_cut']['between_Ma'][1]:.0f} Ma)"))
    h = s2["human"]["eq"]
    top, topw = s2["top_nonhuman_EQ"][0], s2["top_nonhuman_EQ_n3plus"][0]
    rows += [("human EQ", f"{h:.2f}"),
             ("top non-human EQ, single record", f"{top['eq']:.2f} ({top['species']}; {top['ratio_to_human']:.2f} of human)"),
             ("top non-human EQ, species with >= 3 brains", f"{topw['eq']:.2f} ({topw['species']}; {topw['ratio_to_human']:.2f} of human)"),
             ("  same ratios with a 1,350 g human brain", f"{s2['human_brain_1350g_sensitivity']['top_n3plus_ratio']:.2f}-"
                                                         f"{s2['human_brain_1350g_sensitivity']['top_single_record_ratio']:.2f}")]
    for r in s2["lineage_counts"]:
        if r["metric"] == "eq" and r["within_factor"] in (1.5, 2.0, 3.0):
            rows.append((f"species within x{r['within_factor']:g} of human EQ: all records / >= 3 brains",
                         f"{r['n_other_species']} / {r['n_species_n3plus']} (families {r['n_other_families']} / {r['n_families_n3plus']}; "
                         f"lineages separate for 30, 20, 10 Myr: {r['deep30Ma_clades_robust']}, {r['deep20Ma_clades_robust']}, "
                         f"{r['deep10Ma_clades_robust']}; same from species with >= 3 brains only: {r['deep30Ma_clades_n3plus']}, "
                         f"{r['deep20Ma_clades_n3plus']}, {r['deep10Ma_clades_n3plus']})"))
        if r["metric"] == "ncx" and r["within_factor"] in (2.0, 3.0):
            rows.append((f"species within x{r['within_factor']:g} of human cortical neurons (IF; measured where available)",
                         f"{r['n_other_species']} ({r['n_species_measured_count']} with a measured count; families {r['n_other_families']})"))
    es = s2["envelope_slope"]["eq_max"]
    ss = s2["envelope_slope_sensitivity"]
    rows.append(("slope of non-hominin epoch-maximum EQ", f"{es['slope_ln_per_Myr']:.3f} ln units/Myr (doubling {es['doubling_time_Myr']:.0f} Myr); "
                 f"rarefied to n = 10 per epoch {ss['rarefied_n10_per_epoch']['mean']:.3f} "
                 f"({ss['rarefied_n10_per_epoch']['p2_5']:.3f}-{ss['rarefied_n10_per_epoch']['p97_5']:.3f})"))
    R = s2["bm_rates"]["alvarez2022"]
    rows.append(("Brownian rate of ln EQ per Myr: plain / robust / noise-corrected / sister pairs",
                 f"{R['Mammalia_eq']['sigma2']:.4f} / {R['Mammalia_eq']['sigma2_robust']:.4f} / "
                 f"{R['Mammalia_eq']['sigma2_noise_corrected']:.4f} / {R['sister_pairs_eq']['sigma2']:.4f} "
                 f"(tip-error sd {R['Mammalia_eq']['tip_error_sd']:.2f}; noise-corrected 95% interval "
                 f"{R['Mammalia_eq']['sigma2_nc_lo95']:.4f}-{R['Mammalia_eq']['sigma2_nc_hi95']:.4f})"))
    ec = R["estimator_check_simulation"]
    rows.append(("  estimator check (8 simulations at the fitted rate and noise): plain / noise-corrected / sister pairs",
                 f"{ec['plain_range'][0]:.4f}-{ec['plain_range'][1]:.4f} / {ec['whole_tree_noise_corrected_range'][0]:.4f}-"
                 f"{ec['whole_tree_noise_corrected_range'][1]:.4f} / {ec['sister_pair_range'][0]:.4f}-{ec['sister_pair_range'][1]:.4f}"))
    wr_ = s2["waiting_time_range"]
    rows.append(("median first-passage times over all model rows", f"{wr_['min_T50_Myr']:.1f}-{wr_['max_T50_Myr']:.0f} Myr "
                 f"(largest 90th percentile {wr_['max_T90_Myr']:.0f} Myr)"))
    hc = s2["hindcast_driftless_diffusion"]
    rows.append(("hindcast P(any of 5,000 lineages reaches half the human EQ in 61 Myr): plain / noise-corrected rate",
                 f"{hc['standard_N5000_half_human']:.4f} / {hc['noise-corrected_N5000_half_human']:.1e}"))
    hv = s2["hardware_vs_technology"]
    rows.append(("technology hazard, 95% upper limit", f"{hv['rate_upper_limit_95_per_lineage_Myr']:.3f} per lineage-Myr "
                 f"({hv['exposure_lineage_Myr_without_technology']} lineage-Myr, hand-coded); no lower limit"))
    return rows


def main():
    t0 = time.time()
    offline = "--offline" in sys.argv
    os.makedirs(os.path.join(HERE, "..", "data"), exist_ok=True)
    os.makedirs(os.path.join(HERE, "..", "figures"), exist_ok=True)
    if not offline:
        import fetch_data
        print("[1/6] fetching raw data")
        ok = fetch_data.main(force="--force" in sys.argv)
        if not ok:
            print("  some downloads failed; derived tables already in data/ are used for those sources")
    import extract_data
    import build_transcribed
    import build_innovations
    import part1_innovations
    import part2_brains
    print("[2/6] extracting derived tables")
    extract_data.main()
    print("[3/6] hand-transcribed tables")
    build_transcribed.main()
    print("[4/6] innovation tables")
    build_innovations.main()
    print("[5/6] Part 1: innovations")
    s1 = part1_innovations.main()
    print("[6/6] Part 2: brains")
    s2 = part2_brains.main()
    rows = headline_table(s1, s2)
    import csv
    with open(os.path.join(HERE, "..", "data", "headline_numbers.csv"), "w", newline="", encoding="utf8") as f:
        w = csv.writer(f)
        w.writerow(["quantity", "value"])
        w.writerows(rows)
    print("\nHeadline numbers (also written to data/headline_numbers.csv)")
    for q, v in rows:
        print(f"  {q}: {v}")
    print(f"done in {time.time() - t0:.0f} s")


if __name__ == "__main__":
    main()