"""Re-run of Belkin (2026), SSRN 10.2139/ssrn.7644478, from public data.

Our analysis for the factchecked.science page `ap-index-mortality`.
Python 3, standard library only. Reads the downloaded files as data and
never executes anything from them. Run with `python3 -I`:

    python3 -I rerun.py <downloads dir> [<output dir for CSV files>]

<downloads dir> holds, as recorded in downloads/MANIFEST.md:
    gfz-kp/Kp_ap_Ap_SN_F107_since_1932.txt          (GFZ Potsdam, CC BY 4.0)
    silso/TableCyclesMiMa.txt                        (WDC-SILSO, CC BY-NC 4.0)
    un-wpp2024/WPP2024_Demographic_Indicators_Medium.csv.gz
    un-wpp2024/WPP2024_DeathsBySingleAgeSex_Medium_1950-2023.csv.gz
    un-wpp2024/WPP2024_PopulationBySingleAgeSex_Medium_1950-2023.csv.gz
    oulu/Phi_Table_2017.txt                          (cosmic-ray modulation potential, Usoskin et al. 2017)
    worldbank/SP.DYN.CDRT.IN_WLD.json                (optional cross-check; skipped if absent)

With an output folder it also writes annual_series.csv, cycle_year_groups.csv, start_year_scan.csv,
null_c_r12_hist.csv and null_summary.csv there (the numbers behind the page's charts and tiles).
GFZ and the Oulu station update their files, so a fresh download can differ slightly in the latest years from the
copies we used (their checksums are in the page's verification log).

What it does
  1. Annual mean Ap (mean of GFZ daily Ap), world crude death rate (CDR),
     life expectancy at birth (LEx) and a directly age-standardized death
     rate (ASDR, computed here from WPP deaths and mid-year population by
     single year of age), 1950-2023.
  2. Numbers the years of each solar cycle from SILSO's minima under three
     conventions and finds the one that reproduces the preprint's
     R = 0.72 (years 1-12), -0.84 (years 1-4) and 0.89 (years 5-12).
  3. Tests the result: the trend in mortality, the COVID-19 years, the
     famine years, detrending and first differences, ASDR and LEx,
     year-by-year correlations, the choice of split, and three nulls that
     keep the data's structure (shuffling whole solar cycles, shifting the
     solar timeline, and synthetic mortality curves with the same trend and
     year-to-year noise but no solar signal).
"""

import csv
import gzip
import itertools
import json
import math
import os
import random
import statistics
import sys

FIRST, LAST = 1950, 2023
YEARS = list(range(FIRST, LAST + 1))
PAPER = {"r_1_12": 0.72, "r_1_4": -0.84, "r_5_12": 0.89}


# ----------------------------------------------------------------- statistics

def pearson(x, y):
    n = len(x)
    mx = sum(x) / n
    my = sum(y) / n
    sxy = sum((a - mx) * (b - my) for a, b in zip(x, y))
    sxx = sum((a - mx) ** 2 for a in x)
    syy = sum((b - my) ** 2 for b in y)
    return sxy / math.sqrt(sxx * syy)


def _betacf(a, b, x):
    # continued fraction for the incomplete beta function (Lentz)
    tiny = 1e-300
    qab, qap, qam = a + b, a + 1.0, a - 1.0
    c, d = 1.0, 1.0 - qab * x / qap
    d = 1.0 / (d if abs(d) > tiny else tiny)
    h = d
    for m in range(1, 400):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d
        d = 1.0 / (d if abs(d) > tiny else tiny)
        c = 1.0 + aa / c
        c = c if abs(c) > tiny else tiny
        h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d
        d = 1.0 / (d if abs(d) > tiny else tiny)
        c = 1.0 + aa / c
        c = c if abs(c) > tiny else tiny
        dl = d * c
        h *= dl
        if abs(dl - 1.0) < 1e-14:
            break
    return h


def betainc(a, b, x):
    if x <= 0.0:
        return 0.0
    if x >= 1.0:
        return 1.0
    lbt = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log(1.0 - x)
    bt = math.exp(lbt)
    if x < (a + 1.0) / (a + b + 2.0):
        return bt * _betacf(a, b, x) / a
    return 1.0 - bt * _betacf(b, a, 1.0 - x) / b


def p_two_sided(r, n):
    """Two-sided p of a Pearson r with n independent points (t test, n-2 df)."""
    df = n - 2
    if df <= 0:
        return float("nan")
    if abs(r) >= 1.0:
        return 0.0
    t2 = r * r * df / (1.0 - r * r)
    return betainc(df / 2.0, 0.5, df / (df + t2))


def lag1(x):
    m = sum(x) / len(x)
    num = sum((x[i] - m) * (x[i + 1] - m) for i in range(len(x) - 1))
    den = sum((v - m) ** 2 for v in x)
    return num / den


def n_eff(x, y):
    """Effective number of independent pairs for two AR(1)-like series,
    n (1 - r1 r2) / (1 + r1 r2), the usual lag-1 approximation
    (Bretherton et al. 1999; Pyper and Peterman 1998)."""
    r1, r2 = lag1(x), lag1(y)
    n = len(x)
    return n * (1.0 - r1 * r2) / (1.0 + r1 * r2), r1, r2


def solve(a, b):
    """Gaussian elimination with partial pivoting (small dense systems)."""
    n = len(b)
    m = [row[:] + [b[i]] for i, row in enumerate(a)]
    for col in range(n):
        piv = max(range(col, n), key=lambda r: abs(m[r][col]))
        m[col], m[piv] = m[piv], m[col]
        for r in range(col + 1, n):
            f = m[r][col] / m[col][col]
            for c in range(col, n + 1):
                m[r][c] -= f * m[col][c]
    x = [0.0] * n
    for r in range(n - 1, -1, -1):
        x[r] = (m[r][n] - sum(m[r][c] * x[c] for c in range(r + 1, n))) / m[r][r]
    return x


def polyfit_trend(years, values, degree):
    """Least-squares polynomial in scaled time; returns the fitted values."""
    t0 = (years[0] + years[-1]) / 2.0
    h = (years[-1] - years[0]) / 2.0
    s = [(y - t0) / h for y in years]
    k = degree + 1
    ata = [[sum(si ** (i + j) for si in s) for j in range(k)] for i in range(k)]
    atb = [sum(v * si ** i for si, v in zip(s, values)) for i in range(k)]
    coef = solve(ata, atb)
    return [sum(coef[i] * si ** i for i in range(k)) for si in s]


def hp_trend(values, lam):
    """Hodrick-Prescott trend: solve (I + lam * D'D) tau = y, D = second differences."""
    n = len(values)
    a = [[0.0] * n for _ in range(n)]
    for i in range(n):
        a[i][i] = 1.0
    for k in range(n - 2):
        idx = (k, k + 1, k + 2)
        coef = (1.0, -2.0, 1.0)
        for i, ci in zip(idx, coef):
            for j, cj in zip(idx, coef):
                a[i][j] += lam * ci * cj
    return solve(a, list(values))


def fmt_p(p):
    if p != p:
        return "n/a"
    if p < 0.0001:
        return "%.1e" % p
    return "%.4f" % p


# ----------------------------------------------------------------- data

def load_ap(path):
    vals, flags = {}, {}
    with open(path, encoding="utf-8") as f:
        for line in f:
            if line.startswith("#"):
                continue
            p = line.split()
            y = int(p[0])
            ap = int(p[23])
            d = int(p[27])
            if ap < 0:
                continue
            vals.setdefault(y, []).append(ap)
            flags.setdefault(y, []).append(d)
    return vals, flags


def load_minima(path):
    mins = []
    with open(path, encoding="utf-8") as f:
        for line in f:
            p = line.split()
            if len(p) >= 4 and p[0].isdigit() and len(p[0]) == 2:
                mins.append((int(p[0]), int(p[1]), int(p[2])))
    return mins


def load_wpp_indicators(path):
    out = {}
    china = {}
    with gzip.open(path, "rt", encoding="utf-8-sig", newline="") as f:
        for row in csv.DictReader(f):
            t = int(row["Time"])
            if not (FIRST <= t <= LAST):
                continue
            if row["Location"] == "World" and row["LocTypeName"] == "World":
                out[t] = {k: float(row[k]) for k in ("CDR", "LEx", "Deaths", "TPopulation1July", "MedianAgePop")}
            if row["Location"] == "China" and row["LocTypeName"] == "Country/Area":
                china[t] = {k: float(row[k]) for k in ("CDR", "Deaths")}
    return out, china


def load_world_by_age(path, col):
    out = {}
    with gzip.open(path, "rt", encoding="utf-8-sig", newline="") as f:
        for row in csv.DictReader(f):
            if row["Location"] != "World" or row["LocTypeName"] != "World":
                continue
            t = int(row["Time"])
            if FIRST <= t <= LAST:
                out.setdefault(t, {})[int(row["AgeGrpStart"])] = float(row[col])
    return out


def load_worldbank(path):
    with open(path, encoding="utf-8") as f:
        d = json.load(f)
    return {int(x["date"]): x["value"] for x in d[1] if x["value"] is not None}


# ----------------------------------------------------------------- cycle years

def cycle_labels(minima, convention):
    """year -> (cycle, year number). Conventions:
    'min'   year 1 = the calendar year containing SILSO's minimum;
    'after' year 1 = the calendar year after the minimum year (the minimum
            year is the last year of the old cycle);
    'round' year 1 = the minimum year if the minimum falls in January-June,
            otherwise the next year."""
    starts = []
    for c, y, m in minima:
        if convention == "min":
            s = y
        elif convention == "after":
            s = y + 1
        else:
            s = y if m <= 6 else y + 1
        starts.append((c, s))
    lab = {}
    for i, (c, s) in enumerate(starts):
        e = starts[i + 1][1] if i + 1 < len(starts) else 10 ** 4
        for yy in range(s, e):
            lab[yy] = (c, yy - s + 1)
    return lab


def grouped(series, lab, years=None, max_k=12):
    years = YEARS if years is None else years
    g = {}
    for y in years:
        if y in series:
            g.setdefault(lab[y][1], []).append(y)
    ks = [k for k in sorted(g) if k <= max_k]
    means = {k: sum(series[y] for y in g[k]) / len(g[k]) for k in ks}
    return ks, g, means


def three_r(ap, mort, lab, years=None):
    ks, g, apm = grouped(ap, lab, years)
    _, _, mm = grouped(mort, lab, years)
    r = {}
    for name, lo, hi in (("r_1_12", 1, 12), ("r_1_4", 1, 4), ("r_5_12", 5, 12)):
        sel = [k for k in ks if lo <= k <= hi]
        r[name] = (pearson([apm[k] for k in sel], [mm[k] for k in sel]), len(sel))
    return r, ks, g, apm, mm


def show_three(title, r):
    print("  %-58s r(1-12) = %+.3f [n=%d, p=%s]  r(1-4) = %+.3f [n=%d, p=%s]  r(5-12) = %+.3f [n=%d, p=%s]" % (
        title,
        r["r_1_12"][0], r["r_1_12"][1], fmt_p(p_two_sided(*r["r_1_12"])),
        r["r_1_4"][0], r["r_1_4"][1], fmt_p(p_two_sided(*r["r_1_4"])),
        r["r_5_12"][0], r["r_5_12"][1], fmt_p(p_two_sided(*r["r_5_12"]))))


# ----------------------------------------------------------------- main

def main():
    ddir = sys.argv[1]
    outdir = sys.argv[2] if len(sys.argv) > 2 else None
    sys.stdout.reconfigure(encoding="utf-8")

    apd, apflags = load_ap(os.path.join(ddir, "gfz-kp", "Kp_ap_Ap_SN_F107_since_1932.txt"))
    minima = load_minima(os.path.join(ddir, "silso", "TableCyclesMiMa.txt"))
    wpp, china = load_wpp_indicators(os.path.join(ddir, "un-wpp2024", "WPP2024_Demographic_Indicators_Medium.csv.gz"))
    deaths = load_world_by_age(os.path.join(ddir, "un-wpp2024", "WPP2024_DeathsBySingleAgeSex_Medium_1950-2023.csv.gz"), "DeathTotal")
    pop = load_world_by_age(os.path.join(ddir, "un-wpp2024", "WPP2024_PopulationBySingleAgeSex_Medium_1950-2023.csv.gz"), "PopTotal")
    wb_path = os.path.join(ddir, "worldbank", "SP.DYN.CDRT.IN_WLD.json")
    wb = load_worldbank(wb_path) if os.path.exists(wb_path) else None  # optional cross-check

    ap = {y: sum(apd[y]) / len(apd[y]) for y in YEARS}
    cdr = {y: wpp[y]["CDR"] for y in YEARS}
    lex = {y: wpp[y]["LEx"] for y in YEARS}
    neg_lex = {y: -lex[y] for y in YEARS}

    print("=" * 100)
    print("1. DATA")
    print("=" * 100)
    print("Ap: GFZ daily Ap, annual means %d-%d; days per year %d-%d; definitive Kp (flag>=1) on all days: %s" % (
        FIRST, LAST, min(len(apd[y]) for y in YEARS), max(len(apd[y]) for y in YEARS),
        all(min(apflags[y]) >= 1 for y in YEARS)))
    print("SILSO minima (cycle, year, month):", [m for m in minima if m[0] >= 17])

    # checks of the WPP files against each other and against the World Bank
    max_cdr_diff = 0.0
    for y in YEARS:
        d = sum(deaths[y].values())
        p = sum(pop[y].values())
        max_cdr_diff = max(max_cdr_diff, abs(1000.0 * d / p - cdr[y]))
    print("WPP check: deaths/population by single age reproduce WPP's CDR to within %.3f per 1,000 (all years)" % max_cdr_diff)
    print("WPP check: ages in the by-age files: deaths %d-%d, population %d-%d" % (
        min(deaths[2000]), max(deaths[2000]), min(pop[2000]), max(pop[2000])))
    if wb:
        wb_diff = [(y, wb[y] - cdr[y]) for y in YEARS if y in wb]
        print("World Bank SP.DYN.CDRT.IN (World) minus WPP 2024 CDR, %d-%d: max |diff| = %.3f per 1,000 (n=%d years)" % (
            wb_diff[0][0], wb_diff[-1][0], max(abs(d) for _, d in wb_diff), len(wb_diff)))
        big = sorted(wb_diff, key=lambda t: -abs(t[1]))[:6]
        print("  largest differences (year, WB - WPP):", ", ".join("%d %+.2f" % t for t in sorted(big)))
        print("  years with |diff| < 0.1: %d of %d" % (sum(1 for _, d in wb_diff if abs(d) < 0.1), len(wb_diff)))
    else:
        print("World Bank cross-check skipped (worldbank/SP.DYN.CDRT.IN_WLD.json not found; it is optional)")

    # age-standardized death rate (direct method), two standards
    def asdr(std_year):
        w_tot = sum(pop[std_year].values())
        w = {a: v / w_tot for a, v in pop[std_year].items()}
        return {y: 1000.0 * sum(w[a] * deaths[y][a] / pop[y][a] for a in w) for y in YEARS}
    asdr2000 = asdr(2000)
    asdr2023 = asdr(2023)

    print()
    print("Annual series (CDR, ASDR per 1,000; LEx years; China CDR for the famine check)")
    print("year  Ap      CDR     ASDR(2000 std)  ASDR(2023 std)  LEx     medianAge  ChinaCDR")
    for y in YEARS:
        print("%d  %6.2f  %6.3f  %8.3f        %8.3f        %6.2f  %6.2f     %6.2f" % (
            y, ap[y], cdr[y], asdr2000[y], asdr2023[y], lex[y], wpp[y]["MedianAgePop"], china[y]["CDR"]))

    # ------------------------------------------------------------- 2. reproduce
    print()
    print("=" * 100)
    print("2. REPRODUCING THE PREPRINT'S THREE CORRELATIONS (paper: 0.72, -0.84, 0.89)")
    print("=" * 100)
    labs = {c: cycle_labels(minima, c) for c in ("min", "after", "round")}
    for conv in ("min", "after", "round"):
        r, ks, g, apm, mm = three_r(ap, cdr, labs[conv])
        show_three("numbering '%s'" % conv, r)
    lab = labs["after"]
    r_obs, ks, g, apm, mm = three_r(ap, cdr, lab)
    print()
    print("The preprint's numbers are reproduced by numbering 'after' (year 1 = the year after the minimum year).")
    print("Cycle blocks in the record under this numbering:")
    blocks = []
    cur = None
    for y in YEARS:
        c = lab[y][0]
        if cur is None or cur[0] != c:
            cur = [c, []]
            blocks.append(cur)
        cur[1].append(y)
    for c, ys in blocks:
        print("  cycle %d: %d-%d, years %d-%d of the cycle (%d years)" % (c, ys[0], ys[-1], lab[ys[0]][1], lab[ys[-1]][1], len(ys)))
    print()
    print("The 12 data points (group means) behind the correlations:")
    print("  yr  n  mean Ap  mean CDR  mean ASDR2000  mean LEx  mean calendar year  member years")
    _, _, asdrm = grouped(asdr2000, lab)
    _, _, lexm = grouped(lex, lab)
    _, _, yrm = grouped({y: float(y) for y in YEARS}, lab)
    for k in ks:
        print("  %2d %2d  %6.2f   %6.3f    %6.3f         %6.2f    %7.1f            %s" % (
            k, len(g[k]), apm[k], mm[k], asdrm[k], lexm[k], yrm[k], " ".join(str(y) for y in g[k])))
    for name in ("r_1_12", "r_1_4", "r_5_12"):
        rv, n = r_obs[name]
        print("  %s: reproduced %+.4f (paper %+.2f), n = %d points, naive two-sided p = %s" % (
            name, rv, PAPER[name], n, fmt_p(p_two_sided(rv, n))))

    # ------------------------------------------------------------- 3. what drives it
    print()
    print("=" * 100)
    print("3. WHAT DRIVES THE GROUPED CORRELATIONS")
    print("=" * 100)
    show_three("mean calendar year of each group (negated) vs mean Ap", three_r(ap, {y: -float(y) for y in YEARS}, lab)[0])
    for deg in (1, 2, 3, 4):
        tr = polyfit_trend(YEARS, [cdr[y] for y in YEARS], deg)
        trend = dict(zip(YEARS, tr))
        show_three("smooth CDR trend only (polynomial, degree %d)" % deg, three_r(ap, trend, lab)[0])
    # trend fitted without the COVID years, so it carries no 2020-2022 bump
    yrs_nc = [y for y in YEARS if y not in (2020, 2021, 2022)]
    for deg in (3, 4):
        tr = polyfit_trend(yrs_nc, [cdr[y] for y in yrs_nc], deg)
        # fitted on the non-COVID years; the three COVID years are interpolated
        # linearly between the fitted 2019 and 2023 values
        trend = dict(zip(yrs_nc, tr))
        for y in (2020, 2021, 2022):
            trend[y] = trend[2019] + (trend[2023] - trend[2019]) * (y - 2019) / 4.0
        show_three("smooth CDR trend, no COVID bump (degree %d)" % deg, three_r(ap, trend, lab)[0])
    show_three("a straight line falling 0.1 per 1,000 a year", three_r(ap, {y: 20.0 - 0.1 * (y - FIRST) for y in YEARS}, lab)[0])
    # the same grouping applied to series that cannot depend on Ap
    pop_total = {y: wpp[y]["TPopulation1July"] for y in YEARS}
    show_three("world population (negated), no link to Ap possible", three_r(ap, {y: -pop_total[y] for y in YEARS}, lab)[0])
    show_three("world median age (negated)", three_r(ap, {y: -wpp[y]["MedianAgePop"] for y in YEARS}, lab)[0])
    # residual of each group mean from the trend-only group mean
    tr4 = dict(zip(YEARS, polyfit_trend(YEARS, [cdr[y] for y in YEARS], 4)))
    _, _, trm = grouped(tr4, lab)
    print("  group mean CDR minus group mean of the degree-4 trend (what the trend does not explain):")
    print("    " + "  ".join("yr%d %+.3f" % (k, mm[k] - trm[k]) for k in ks))
    hp100 = dict(zip(YEARS, hp_trend([cdr[y] for y in YEARS], 100.0)))
    print("  CDR residuals, degree-4 polynomial / HP(100), early and late years:")
    print("    " + "  ".join("%d %+.2f/%+.2f" % (y, cdr[y] - tr4[y], cdr[y] - hp100[y]) for y in list(range(1950, 1966)) + list(range(2017, 2024))))
    print("  groups 1-10 only (drops the thin groups 11 and 12, 3 and 2 years):")
    sel = [k for k in ks if k <= 10]
    r10 = pearson([apm[k] for k in sel], [mm[k] for k in sel])
    sel2 = [k for k in ks if 5 <= k <= 10]
    r510 = pearson([apm[k] for k in sel2], [mm[k] for k in sel2])
    print("    r(1-10) = %+.3f (n=%d, p=%s); r(5-10) = %+.3f (n=%d, p=%s)" % (
        r10, len(sel), fmt_p(p_two_sided(r10, len(sel))), r510, len(sel2), fmt_p(p_two_sided(r510, len(sel2)))))
    print("  start and end years of the record (CDR, numbering 'after'):")
    for a0, a1 in ((1950, 2023), (1951, 2023), (1952, 2023), (1953, 2023), (1954, 2023), (1955, 2023), (1956, 2023),
                   (1958, 2023), (1960, 2023), (1962, 2023), (1965, 2023), (1970, 2023), (1975, 2023), (1980, 2023),
                   (1950, 2019), (1950, 2015), (1950, 2010), (1950, 2000),
                   (1950, 1986), (1987, 2023)):
        rr_ = three_r(ap, cdr, lab, list(range(a0, a1 + 1)))[0]
        print("    %d-%d: r(1-12) = %+.3f  r(1-4) = %+.3f  r(5-12) = %+.3f" % (a0, a1, rr_["r_1_12"][0], rr_["r_1_4"][0], rr_["r_5_12"][0]))
    for label_, spans in (("every start year 1950-1980, ending 2023", [(a, 2023) for a in range(1950, 1981)]),
                          ("every end year 1993-2023, starting 1950", [(1950, b) for b in range(1993, 2024)])):
        rs_ = [three_r(ap, cdr, lab, list(range(a, b + 1)))[0] for a, b in spans]
        v12 = [x["r_1_12"][0] for x in rs_]
        v14 = [x["r_1_4"][0] for x in rs_]
        v512 = [x["r_5_12"][0] for x in rs_]
        print("    %s (%d spans): r(1-12) %+.2f to %+.2f; r(1-4) %+.2f to %+.2f (positive in %d); r(5-12) %+.2f to %+.2f" % (
            label_, len(spans), min(v12), max(v12), min(v14), max(v14), sum(1 for v in v14 if v > 0), min(v512), max(v512)))
        if label_.startswith("every start"):
            line_ = {y: 20.0 - 0.1 * (y - FIRST) for y in YEARS}
            scan_rows = [(a, b, x["r_1_12"][0], x["r_1_4"][0], x["r_5_12"][0],
                          three_r(ap, line_, lab, list(range(a, b + 1)))[0]["r_1_12"][0]) for (a, b), x in zip(spans, rs_)]
            print("      the same scan for a straight falling line: r(1-12) %+.2f to %+.2f; across the %d start years the two "
                  "curves of r(1-12) correlate at %.3f" % (min(r[5] for r in scan_rows), max(r[5] for r in scan_rows), len(scan_rows),
                                                          pearson([r[2] for r in scan_rows], [r[5] for r in scan_rows])))

    # ------------------------------------------------------------- 4. robustness
    print()
    print("=" * 100)
    print("4. ROBUSTNESS: SAME METHOD, OTHER CHOICES")
    print("=" * 100)
    show_three("as published: CDR 1950-2023", r_obs)
    show_three("without the COVID-19 years 2020-2022", three_r(ap, cdr, lab, [y for y in YEARS if y not in (2020, 2021, 2022)])[0])
    show_three("without 2020-2023", three_r(ap, cdr, lab, [y for y in YEARS if y < 2020])[0])
    show_three("without the famine years 1959-1961", three_r(ap, cdr, lab, [y for y in YEARS if y not in (1959, 1960, 1961)])[0])
    show_three("without 1959-1961 and 2020-2022", three_r(ap, cdr, lab, [y for y in YEARS if y not in (1959, 1960, 1961, 2020, 2021, 2022)])[0])
    show_three("complete cycles only, 1955-2019 (cycles 19-24)", three_r(ap, cdr, lab, list(range(1955, 2020)))[0])
    show_three("1960-2023 (the World Bank's span)", three_r(ap, cdr, lab, list(range(1960, 2024)))[0])
    show_three("1970-2023", three_r(ap, cdr, lab, list(range(1970, 2024)))[0])
    show_three("age-standardized death rate (2000 standard)", three_r(ap, asdr2000, lab)[0])
    show_three("age-standardized death rate (2023 standard)", three_r(ap, asdr2023, lab)[0])
    show_three("life expectancy at birth, negated", three_r(ap, neg_lex, lab)[0])
    shocks = (1959, 1960, 1961, 2020, 2021, 2022)
    for name, series in (("CDR", cdr), ("ASDR2000", asdr2000), ("-LEx", neg_lex)):
        vals = [series[y] for y in YEARS]
        for deg in (2, 3, 4):
            tr = polyfit_trend(YEARS, vals, deg)
            res = {y: v - t for y, v, t in zip(YEARS, vals, tr)}
            show_three("%s minus polynomial trend (degree %d)" % (name, deg), three_r(ap, res, lab)[0])
        for lam in (100.0, 400.0):
            tr = hp_trend(vals, lam)
            res = {y: v - t for y, v, t in zip(YEARS, vals, tr)}
            show_three("%s minus Hodrick-Prescott trend (lambda %d)" % (name, lam), three_r(ap, res, lab)[0])
        # trend fitted without the shock years, residuals scored without them
        yrs_ns = [y for y in YEARS if y not in shocks]
        tr = polyfit_trend(yrs_ns, [series[y] for y in yrs_ns], 4)
        res = {y: series[y] - t for y, t in zip(yrs_ns, tr)}
        show_three("%s minus degree-4 trend, without 1959-61 and 2020-22" % name, three_r(ap, res, lab, yrs_ns)[0])
        if name == "CDR":
            _, g_ns, rgm = grouped(res, lab, yrs_ns)
            mean_rate = statistics.mean(series[y] for y in yrs_ns)
            print("    its 12 residual group means span %+.3f to %+.3f per 1,000, %.1f%% of the average rate (%.2f per 1,000)" % (
                min(rgm.values()), max(rgm.values()), 100.0 * (max(rgm.values()) - min(rgm.values())) / mean_rate, mean_rate))
            print("    residual group means: " + "  ".join("yr%d %+.3f" % (k, rgm[k]) for k in sorted(rgm)))
            big1 = max(g_ns[1], key=lambda y: abs(res[y]))
            print("    the year-1 group's largest departure from the trend: %d (%+.2f per 1,000)" % (big1, res[big1]))
            yrs_ns1 = [y for y in yrs_ns if y != big1]
            tr1 = polyfit_trend(yrs_ns1, [series[y] for y in yrs_ns1], 4)
            res1 = {y: series[y] - t for y, t in zip(yrs_ns1, tr1)}
            show_three("CDR minus degree-4 trend, without the shocks and %d" % big1, three_r(ap, res1, lab, yrs_ns1)[0])
            rgm_cdr, mean_rate_ns = dict(rgm), mean_rate  # for the scale check
        diff = {y: series[y] - series[y - 1] for y in YEARS[1:]}
        dap = {y: ap[y] - ap[y - 1] for y in YEARS[1:]}
        # Ap levels against changes in mortality: kept for the record; it is not a test of a link between levels
        show_three("%s changes against Ap levels (1951-2023; not a level test)" % name, three_r(ap, diff, lab, YEARS[1:])[0])
        yrs_d = [y for y in YEARS[1:] if y not in shocks and y - 1 not in shocks]
        show_three("%s changes against Ap levels, without changes into or out of the shocks" % name, three_r(ap, diff, lab, yrs_d)[0])
        # like for like: changes in mortality against changes in Ap
        show_three("%s changes against changes in Ap (1951-2023)" % name, three_r(dap, diff, lab, YEARS[1:])[0])
        logr = {y: math.log(abs(series[y]) / abs(series[y - 1])) for y in YEARS[1:]}
        show_three("%s year-on-year log change against Ap levels (1951-2023)" % name, three_r(ap, logr, lab, YEARS[1:])[0])
    print()
    print("Other year numberings (CDR 1950-2023):")
    for conv in ("min", "round"):
        show_three("numbering '%s'" % conv, three_r(ap, cdr, labs[conv])[0])
    print()
    print("Where to split the cycle (paper: growth = years 1-4, decline = 5-12). CDR, numbering 'after':")
    for k in range(2, 11):
        lo = [kk for kk in ks if kk <= k]
        hi = [kk for kk in ks if kk > k]
        r_lo = pearson([apm[x] for x in lo], [mm[x] for x in lo]) if len(lo) >= 3 else float("nan")
        r_hi = pearson([apm[x] for x in hi], [mm[x] for x in hi]) if len(hi) >= 3 else float("nan")
        print("  years 1-%d: r = %+.3f (n=%d)   years %d-12: r = %+.3f (n=%d)" % (k, r_lo, len(lo), k + 1, r_hi, len(hi)))

    # ------------------------------------------------------------- 5. year by year
    print()
    print("=" * 100)
    print("5. YEAR BY YEAR (no grouping)")
    print("=" * 100)
    xs = [ap[y] for y in YEARS]
    for name, series in (("CDR", cdr), ("ASDR2000", asdr2000), ("-LEx", neg_lex)):
        ys = [series[y] for y in YEARS]
        rr = pearson(xs, ys)
        ne, r1x, r1y = n_eff(xs, ys)
        print("  Ap vs %-9s levels: r = %+.3f, n = %d, naive p = %s; lag-1 autocorrelations %.2f and %.2f, effective n = %.1f, p = %s" % (
            name, rr, len(xs), fmt_p(p_two_sided(rr, len(xs))), r1x, r1y, ne, fmt_p(p_two_sided(rr, max(3, int(ne))))))
        dx = [ap[y] - ap[y - 1] for y in YEARS[1:]]
        dy = [series[y] - series[y - 1] for y in YEARS[1:]]
        rd = pearson(dx, dy)
        print("  dAp vs d%-8s first differences: r = %+.3f, n = %d, p = %s" % (name, rd, len(dx), fmt_p(p_two_sided(rd, len(dx)))))
        ra = pearson([ap[y] for y in YEARS[1:]], dy)
        print("  Ap vs d%-9s (level of Ap vs change in mortality): r = %+.3f, n = %d, p = %s" % (name, ra, len(dy), fmt_p(p_two_sided(ra, len(dy)))))
        for deg in (2, 3, 4):
            tr = polyfit_trend(YEARS, ys, deg)
            res = [v - t for v, t in zip(ys, tr)]
            rres = pearson(xs, res)
            ne, r1x, r1y = n_eff(xs, res)
            print("  Ap vs %-9s minus degree-%d trend: r = %+.3f, n = %d, naive p = %s, effective n = %.1f, p = %s" % (
                name, deg, rres, len(xs), fmt_p(p_two_sided(rres, len(xs))), ne, fmt_p(p_two_sided(rres, max(3, int(ne))))))
        tr = hp_trend(ys, 100.0)
        res = [v - t for v, t in zip(ys, tr)]
        rres = pearson(xs, res)
        ne, r1x, r1y = n_eff(xs, res)
        print("  Ap vs %-9s minus HP trend (lambda 100): r = %+.3f, n = %d, naive p = %s, effective n = %.1f, p = %s" % (
            name, rres, len(xs), fmt_p(p_two_sided(rres, len(xs))), ne, fmt_p(p_two_sided(rres, max(3, int(ne))))))
        sel = [i for i, y in enumerate(YEARS) if y not in shocks]
        rres2 = pearson([xs[i] for i in sel], [res[i] for i in sel])
        print("  Ap vs %-9s minus HP trend (lambda 100), without 1959-61 and 2020-22: r = %+.3f, n = %d, naive p = %s" % (
            name, rres2, len(sel), fmt_p(p_two_sided(rres2, len(sel)))))
        sel = [y for y in YEARS[1:] if y not in shocks and y - 1 not in shocks]
        rd2 = pearson([ap[y] - ap[y - 1] for y in sel], [series[y] - series[y - 1] for y in sel])
        print("  dAp vs d%-8s first differences, without changes into or out of 1959-61 and 2020-22: r = %+.3f, n = %d, p = %s" % (
            name, rd2, len(sel), fmt_p(p_two_sided(rd2, len(sel)))))
        # growth phase only (paper's years 1-4) and years 5-12: Ap against detrended mortality, year by year
        hp_res = dict(zip(YEARS, [v - t for v, t in zip(ys, hp_trend(ys, 100.0))]))
        for lo, hi in ((1, 4), (5, 12)):
            sel = [y for y in YEARS if lo <= lab[y][1] <= hi]
            rr2 = pearson([ap[y] for y in sel], [hp_res[y] for y in sel])
            sel_ns = [y for y in sel if y not in shocks]
            rr3 = pearson([ap[y] for y in sel_ns], [hp_res[y] for y in sel_ns])
            print("  Ap vs %s minus HP trend, years %d-%d only: r = %+.3f, n = %d, p = %s; without the shock years r = %+.3f, n = %d, p = %s" % (
                name, lo, hi, rr2, len(sel), fmt_p(p_two_sided(rr2, len(sel))), rr3, len(sel_ns), fmt_p(p_two_sided(rr3, len(sel_ns)))))
        # growth phase only (paper's years 1-4) and years 5-12, first differences
        for lo, hi in ((1, 4), (5, 12)):
            sel = [y for y in YEARS[1:] if lo <= lab[y][1] <= hi]
            rr2 = pearson([ap[y] - ap[y - 1] for y in sel], [series[y] - series[y - 1] for y in sel])
            sel_ns = [y for y in sel if y not in shocks and y - 1 not in shocks]
            rr3 = pearson([ap[y] - ap[y - 1] for y in sel_ns], [series[y] - series[y - 1] for y in sel_ns])
            print("  first differences, years %d-%d only: r = %+.3f, n = %d, p = %s; without the shock years r = %+.3f, n = %d, p = %s" % (
                lo, hi, rr2, len(sel), fmt_p(p_two_sided(rr2, len(sel))), rr3, len(sel_ns), fmt_p(p_two_sided(rr3, len(sel_ns)))))
    # 5b. the same phase-specific tests with a trend that the shocks cannot pull: the famine and
    #     COVID years are replaced by straight-line interpolation before the trend is fitted, and
    #     left out of the correlations; "within cycle" also removes each cycle's own mean
    print()
    print("  5b. Phase by phase, trend fitted with 1959-61 and 2020-22 interpolated (shock years left out of the correlations):")
    for name, series in (("CDR", cdr), ("ASDR2000", asdr2000), ("-LEx", neg_lex)):
        sf = dict(series)
        for a0, b0, ys_ in ((1958, 1962, (1959, 1960, 1961)), (2019, 2023, (2020, 2021, 2022))):
            for y in ys_:
                sf[y] = series[a0] + (series[b0] - series[a0]) * (y - a0) / (b0 - a0)
        vals_sf = [sf[y] for y in YEARS]
        for kind, tr in (("HP 100", hp_trend(vals_sf, 100.0)), ("HP 400", hp_trend(vals_sf, 400.0)),
                         ("degree 3", polyfit_trend(YEARS, vals_sf, 3)), ("degree 4", polyfit_trend(YEARS, vals_sf, 4))):
            res = {y: series[y] - t for y, t in zip(YEARS, tr)}
            parts = []
            for lo, hi in ((1, 4), (5, 12), (1, 12)):
                sel = [y for y in YEARS if lo <= lab[y][1] <= hi and y not in shocks]
                rr4 = pearson([ap[y] for y in sel], [res[y] for y in sel])
                byc = {}
                for y in sel:
                    byc.setdefault(lab[y][0], []).append(y)
                xa, xr = [], []
                for c, ys_ in byc.items():
                    if len(ys_) < 2:
                        continue
                    ma = statistics.mean(ap[y] for y in ys_)
                    mr = statistics.mean(res[y] for y in ys_)
                    xa.extend(ap[y] - ma for y in ys_)
                    xr.extend(res[y] - mr for y in ys_)
                rw = pearson(xa, xr)
                parts.append("years %d-%d r = %+.3f (n=%d, p=%s), within cycle %+.3f" % (lo, hi, rr4, len(sel), fmt_p(p_two_sided(rr4, len(sel))), rw))
            print("    %-8s %-8s %s" % (name, kind, " | ".join(parts)))

    # solar-minimum years (cosmic rays peak near minimum) against the others, detrended ASDR
    res = dict(zip(YEARS, [v - t for v, t in zip([asdr2000[y] for y in YEARS], polyfit_trend(YEARS, [asdr2000[y] for y in YEARS], 4))]))
    minyears = [y for (c, y, m) in minima if FIRST <= y <= LAST]
    near = set()
    for y in minyears:
        near.update((y - 1, y, y + 1))
    near &= set(YEARS)
    a = [res[y] for y in YEARS if y in near]
    b = [res[y] for y in YEARS if y not in near]
    print()
    print("  Solar minimum years %s and the years either side (n=%d) vs all other years (n=%d):" % (minyears, len(a), len(b)))
    print("    ASDR (2000 standard) minus degree-4 trend: mean %+.3f vs %+.3f per 1,000" % (statistics.mean(a), statistics.mean(b)))
    sn = {}
    with open(os.path.join(ddir, "gfz-kp", "Kp_ap_Ap_SN_F107_since_1932.txt"), encoding="utf-8") as f:
        for line in f:
            if line.startswith("#"):
                continue
            p = line.split()
            y = int(p[0])
            s = int(p[24])
            if s >= 0:
                sn.setdefault(y, []).append(s)
    sny = {y: sum(sn[y]) / len(sn[y]) for y in YEARS}
    rs = pearson([sny[y] for y in YEARS], [res[y] for y in YEARS])
    ne, _, _ = n_eff([sny[y] for y in YEARS], [res[y] for y in YEARS])
    print("    annual sunspot number vs detrended ASDR: r = %+.3f, n = %d, naive p = %s, effective n = %.1f, p = %s" % (
        rs, len(YEARS), fmt_p(p_two_sided(rs, len(YEARS))), ne, fmt_p(p_two_sided(rs, max(3, int(ne))))))
    # cosmic-ray modulation potential phi (Usoskin et al. 2017 reference reconstruction, Oulu):
    # low phi = weak solar modulation = more galactic cosmic rays at Earth
    phi = {}
    raw = open(os.path.join(ddir, "oulu", "Phi_Table_2017.txt"), "rb").read()
    txt = raw.decode("utf-16")
    for line in txt.splitlines():
        p = line.split()
        if len(p) >= 14 and p[0].isdigit() and len(p[0]) == 4:
            try:
                phi[int(p[0])] = float(p[13])
            except ValueError:
                pass
    yrs_phi = [y for y in YEARS if y in phi]
    print("    cosmic-ray modulation potential phi (annual, MV), %d-%d: %d years; lowest phi (most cosmic rays) in %s" % (
        yrs_phi[0], yrs_phi[-1], len(yrs_phi), ", ".join(str(y) for y in sorted(yrs_phi, key=lambda y: phi[y])[:5])))
    print("    r(phi, Ap) = %+.3f; r(phi, sunspot number) = %+.3f (n = %d)" % (
        pearson([phi[y] for y in yrs_phi], [ap[y] for y in yrs_phi]), pearson([phi[y] for y in yrs_phi], [sny[y] for y in yrs_phi]), len(yrs_phi)))
    print("    highest phi (fewest cosmic rays) in %s; phi in 1959-61: %s; in 2019-21: %s" % (
        ", ".join("%d (%.0f)" % (y, phi[y]) for y in sorted(yrs_phi, key=lambda y: -phi[y])[:8]),
        ", ".join("%.0f" % phi[y] for y in (1959, 1960, 1961)), ", ".join("%.0f" % phi[y] for y in (2019, 2020, 2021))))
    famine, pandemic = (1959, 1960, 1961), (2020, 2021, 2022)
    variants_phi = {"shock years out": [], "both left in": [], "only the pandemic years left in": [], "only the famine years left in": []}
    for name, series in (("CDR", cdr), ("ASDR2000", asdr2000), ("-LEx", neg_lex)):
        sf = dict(series)
        for a0, b0, ys_ in ((1958, 1962, (1959, 1960, 1961)), (2019, 2023, (2020, 2021, 2022))):
            for y in ys_:
                sf[y] = series[a0] + (series[b0] - series[a0]) * (y - a0) / (b0 - a0)
        for kind, tr in (("HP 100", hp_trend([sf[y] for y in YEARS], 100.0)), ("degree 4", polyfit_trend(YEARS, [sf[y] for y in YEARS], 4))):
            resd = {y: series[y] - t for y, t in zip(YEARS, tr)}
            sel = [y for y in yrs_phi if y not in shocks]
            rp = pearson([phi[y] for y in sel], [resd[y] for y in sel])
            ne, _, _ = n_eff([phi[y] for y in sel], [resd[y] for y in sel])
            sel_all = yrs_phi
            rp_all = pearson([phi[y] for y in sel_all], [resd[y] for y in sel_all])
            print("    phi vs %-8s minus %-8s trend (shock years out): r = %+.3f, n = %d, naive p = %s, effective n = %.1f, p = %s; with shock years r = %+.3f (n = %d)" % (
                name, kind, rp, len(sel), fmt_p(p_two_sided(rp, len(sel))), ne, fmt_p(p_two_sided(rp, max(3, int(ne)))), rp_all, len(sel_all)))
            sel_pan = [y for y in yrs_phi if y not in famine]
            sel_fam = [y for y in yrs_phi if y not in pandemic]
            r_pan = pearson([phi[y] for y in sel_pan], [resd[y] for y in sel_pan])
            r_fam = pearson([phi[y] for y in sel_fam], [resd[y] for y in sel_fam])
            print("      only the pandemic years left in: r = %+.3f (n = %d, p = %s); only the famine years left in: r = %+.3f (n = %d)" % (
                r_pan, len(sel_pan), fmt_p(p_two_sided(r_pan, len(sel_pan))), r_fam, len(sel_fam)))
            variants_phi["shock years out"].append((rp, len(sel)))
            variants_phi["both left in"].append((rp_all, len(sel_all)))
            variants_phi["only the pandemic years left in"].append((r_pan, len(sel_pan)))
            variants_phi["only the famine years left in"].append((r_fam, len(sel_fam)))
    for k, v in variants_phi.items():
        print("    phi against detrended mortality, %s: r %+.2f to %+.2f across the 6 versions (n = %d)" % (
            k, min(x for x, _ in v), max(x for x, _ in v), v[0][1]))

    # ------------------------------------------------------------- 6. nulls
    print()
    print("=" * 100)
    print("6. NULLS THAT KEEP THE DATA'S STRUCTURE (numbering 'after', 1950-2023)")
    print("=" * 100)

    def stats_from_sequence(seq_labels, seq_ap, mort):
        # seq_* are aligned with YEARS
        ssum, scnt, msum = {}, {}, {}
        for i, y in enumerate(YEARS):
            k = seq_labels[i]
            if k > 12:
                continue
            ssum[k] = ssum.get(k, 0.0) + seq_ap[i]
            msum[k] = msum.get(k, 0.0) + mort[y]
            scnt[k] = scnt.get(k, 0) + 1
        kk = sorted(ssum)
        a_ = {k: ssum[k] / scnt[k] for k in kk}
        m_ = {k: msum[k] / scnt[k] for k in kk}
        out = {}
        for name, lo, hi in (("r_1_12", 1, 12), ("r_1_4", 1, 4), ("r_5_12", 5, 12)):
            sel = [k for k in kk if lo <= k <= hi]
            out[name] = pearson([a_[k] for k in sel], [m_[k] for k in sel])
        return out

    def summarize(title, obs, results):
        n = len(results)
        ge12 = sum(1 for r in results if r["r_1_12"] >= obs["r_1_12"] - 1e-12)
        le14 = sum(1 for r in results if r["r_1_4"] <= obs["r_1_4"] + 1e-12)
        ge512 = sum(1 for r in results if r["r_5_12"] >= obs["r_5_12"] - 1e-12)
        joint = sum(1 for r in results if r["r_1_12"] >= obs["r_1_12"] - 1e-12 and r["r_1_4"] <= obs["r_1_4"] + 1e-12 and r["r_5_12"] >= obs["r_5_12"] - 1e-12)
        r12 = sorted(r["r_1_12"] for r in results)
        print("  %s: %d alternatives; observed r(1-12) %+.3f, r(1-4) %+.3f, r(5-12) %+.3f" % (
            title, n, obs["r_1_12"], obs["r_1_4"], obs["r_5_12"]))
        print("    r(1-12) at least as high in %d (p = %.4f); median r(1-12) = %+.3f, 5th-95th percentile %+.3f to %+.3f" % (
            ge12, ge12 / n, statistics.median(r12), r12[int(0.05 * (n - 1))], r12[int(0.95 * (n - 1))]))
        print("    r(1-4) at least as negative in %d (p = %.4f); r(5-12) at least as high in %d (p = %.4f); all three at once in %d (p = %.4f)" % (
            le14, le14 / n, ge512, ge512 / n, joint, joint / n))
        return {"runs": n, "p_r_1_12": ge12 / n, "p_r_1_4": le14 / n, "p_r_5_12": ge512 / n, "p_all_three": joint / n,
                "median_r_1_12": statistics.median(r12), "p05_r_1_12": r12[int(0.05 * (n - 1))], "p95_r_1_12": r12[int(0.95 * (n - 1))]}

    tr4_all = dict(zip(YEARS, polyfit_trend(YEARS, [cdr[y] for y in YEARS], 4)))
    det4 = {y: cdr[y] - tr4_all[y] for y in YEARS}
    det4_ns = {y: (0.0 if y in shocks else det4[y]) for y in YEARS}
    line = {y: 20.0 - 0.1 * (y - FIRST) for y in YEARS}
    test_series = (
        ("CDR as published", cdr),
        ("smooth trend only (degree-4 fit to CDR)", tr4_all),
        ("straight falling line", line),
        ("CDR minus degree-4 trend", det4),
        ("CDR minus degree-4 trend, 1959-61 and 2020-22 set to the trend", det4_ns),
        ("age-standardized death rate (2000 standard)", asdr2000),
    )

    # (a) shuffle whole solar cycles: each cycle keeps its own sequence of year numbers and Ap values;
    #     mortality stays where it is in time
    blk = [[(lab[y][1], ap[y]) for y in ys] for c, ys in blocks]
    print("(a) Shuffle the order of the %d solar-cycle blocks (all %d orderings; the observed one included)" % (
        len(blk), math.factorial(len(blk))))
    for title, series in test_series:
        obs = {k: v[0] for k, v in three_r(ap, series, lab)[0].items()}
        res_all = []
        for order in itertools.permutations(range(len(blk))):
            seq = [item for i in order for item in blk[i]]
            res_all.append(stats_from_sequence([s[0] for s in seq], [s[1] for s in seq], series))
        summarize(title, obs, res_all)
    print("(a') Shuffle only the 6 complete cycles, partial cycles 18 and 25 kept at the ends (720 orderings)")
    inner = list(range(1, len(blk) - 1))
    for title, series in test_series[:2]:
        obs = {k: v[0] for k, v in three_r(ap, series, lab)[0].items()}
        res_in = []
        for order in itertools.permutations(inner):
            o = [0] + list(order) + [len(blk) - 1]
            seq = [item for i in o for item in blk[i]]
            res_in.append(stats_from_sequence([s[0] for s in seq], [s[1] for s in seq], series))
        summarize(title, obs, res_in)

    # (b) shift the solar timeline (year numbers with their Ap values) against mortality, wrapping round
    seq_lab = [lab[y][1] for y in YEARS]
    seq_ap = [ap[y] for y in YEARS]
    print("(b) Shift the solar timeline against mortality by 0-73 years, wrapping round (74 alignments)")
    for title, series in test_series:
        obs = {k: v[0] for k, v in three_r(ap, series, lab)[0].items()}
        res_sh = []
        for k in range(len(YEARS)):
            sl = seq_lab[-k:] + seq_lab[:-k] if k else seq_lab
            sa = seq_ap[-k:] + seq_ap[:-k] if k else seq_ap
            res_sh.append(stats_from_sequence(sl, sa, series))
        summarize(title, obs, res_sh)

    # (c) synthetic mortality with no solar signal: CDR's own smooth trend plus noise like its residuals;
    #     the real Ap and the real year numbering are kept
    print("(c) Synthetic mortality curves with CDR's trend and year-to-year noise but no solar signal (real Ap, real numbering)")
    obs = {k: v[0] for k, v in r_obs.items()}
    rng = random.Random(20261024)
    vals = [cdr[y] for y in YEARS]
    null_r12 = {}  # the r(1-12) values of each set of runs, for the histogram on the page
    null_summary = {}  # p-values and percentiles of each set of runs, for the page
    for deg in (3, 4):
        tr = polyfit_trend(YEARS, vals, deg)
        resid = [v - t for v, t in zip(vals, tr)]
        phi = lag1(resid)
        sd = statistics.pstdev(resid)
        innov = sd * math.sqrt(max(1e-9, 1 - phi * phi))
        sims = []
        for _ in range(20000):
            e = rng.gauss(0.0, sd)
            synth = {}
            for i, y in enumerate(YEARS):
                if i:
                    e = phi * e + rng.gauss(0.0, innov)
                synth[y] = tr[i] + e
            sims.append(stats_from_sequence(seq_lab, seq_ap, synth))
        null_summary["deg%d_ar1" % deg] = summarize("degree-%d trend + AR(1) noise (lag-1 %.2f, sd %.3f), 20,000 runs" % (deg, phi, sd), obs, sims)
        null_r12["deg%d_ar1" % deg] = [r["r_1_12"] for r in sims]
        # the real residuals, resampled in circular blocks of 5 years (keeps spikes like the famine and COVID)
        sims = []
        nY = len(YEARS)
        for _ in range(20000):
            e = []
            while len(e) < nY:
                s0 = rng.randrange(nY)
                e.extend(resid[(s0 + j) % nY] for j in range(5))
            synth = {y: tr[i] + e[i] for i, y in enumerate(YEARS)}
            sims.append(stats_from_sequence(seq_lab, seq_ap, synth))
        null_summary["deg%d_blocks" % deg] = summarize("degree-%d trend + its own residuals in 5-year blocks, 20,000 runs" % deg, obs, sims)
        null_r12["deg%d_blocks" % deg] = [r["r_1_12"] for r in sims]
    print("  across the four versions: r(1-12) p %.2f-%.2f, r(1-4) p %.2f-%.2f, r(5-12) p %.3f-%.3f, all three p %.3f-%.3f" % (
        min(s["p_r_1_12"] for s in null_summary.values()), max(s["p_r_1_12"] for s in null_summary.values()),
        min(s["p_r_1_4"] for s in null_summary.values()), max(s["p_r_1_4"] for s in null_summary.values()),
        min(s["p_r_5_12"] for s in null_summary.values()), max(s["p_r_5_12"] for s in null_summary.values()),
        min(s["p_all_three"] for s in null_summary.values()), max(s["p_all_three"] for s in null_summary.values())))

    # (d) trend-free noise against the detrended statistic: how often noise like the death rate's own departures
    #     from its trend, with no solar signal, gives a grouped correlation as high as the detrended one
    print("(d) Trend-free noise like the death rate's departures from its trend, against the detrended correlations")
    obs_d = {k: v[0] for k, v in three_r(ap, det4, lab)[0].items()}
    resid4 = [det4[y] for y in YEARS]
    phi_d, sd_d = lag1(resid4), statistics.pstdev(resid4)
    innov_d = sd_d * math.sqrt(max(1e-9, 1 - phi_d * phi_d))
    sims = []
    for _ in range(20000):
        e = rng.gauss(0.0, sd_d)
        synth = {}
        for i, y in enumerate(YEARS):
            if i:
                e = phi_d * e + rng.gauss(0.0, innov_d)
            synth[y] = e
        sims.append(stats_from_sequence(seq_lab, seq_ap, synth))
    null_summary["detrended_ar1"] = summarize("CDR minus degree-4 trend; null: AR(1) noise (lag-1 %.2f, sd %.3f), 20,000 runs" % (phi_d, sd_d), obs_d, sims)
    sims = []
    nY = len(YEARS)
    for _ in range(20000):
        e = []
        while len(e) < nY:
            s0 = rng.randrange(nY)
            e.extend(resid4[(s0 + j) % nY] for j in range(5))
        synth = {y: e[i] for i, y in enumerate(YEARS)}
        sims.append(stats_from_sequence(seq_lab, seq_ap, synth))
    null_summary["detrended_blocks"] = summarize("CDR minus degree-4 trend; null: its own departures in 5-year blocks, 20,000 runs", obs_d, sims)
    yrs_ns = [y for y in YEARS if y not in shocks]
    tr_ns = polyfit_trend(yrs_ns, [cdr[y] for y in yrs_ns], 4)
    res_ns = [cdr[y] - t for y, t in zip(yrs_ns, tr_ns)]
    obs_ns = {k: v[0] for k, v in three_r(ap, dict(zip(yrs_ns, res_ns)), lab, yrs_ns)[0].items()}
    phi_n, sd_n = lag1(res_ns), statistics.pstdev(res_ns)
    innov_n = sd_n * math.sqrt(max(1e-9, 1 - phi_n * phi_n))
    sims = []
    for _ in range(20000):
        e = rng.gauss(0.0, sd_n)
        synth = {}
        for i, y in enumerate(yrs_ns):
            if i:
                e = phi_n * e + rng.gauss(0.0, innov_n)
            synth[y] = e
        sims.append({k: v[0] for k, v in three_r(ap, synth, lab, yrs_ns)[0].items()})
    null_summary["detrended_noshock_ar1"] = summarize(
        "CDR minus degree-4 trend fitted without 1959-61 and 2020-22, those years left out; null: AR(1) noise (lag-1 %.2f, sd %.3f), 20,000 runs" % (
            phi_n, sd_n), obs_ns, sims)

    # ------------------------------------------------------------- 7. scale check
    print()
    print("=" * 100)
    print("7. SCALE CHECK AGAINST THE 263-CITY STUDY (Zilli Vieira et al. 2019: +0.13% to +0.47% deaths per SD of daily Kp)")
    print("=" * 100)
    kpd = {}
    with open(os.path.join(ddir, "gfz-kp", "Kp_ap_Ap_SN_F107_since_1932.txt"), encoding="utf-8") as f:
        for line in f:
            if line.startswith("#"):
                continue
            p = line.split()
            y = int(p[0])
            if not (FIRST <= y <= LAST):
                continue
            kps = [float(v) for v in p[7:15]]
            if min(kps) < 0:
                continue
            kpd.setdefault(y, []).append(sum(kps) / 8.0)
    allkp = [v for y in YEARS for v in kpd[y]]
    sdk = statistics.pstdev(allkp)
    kpy = {y: statistics.mean(kpd[y]) for y in YEARS}
    _, _, kpm = grouped(kpy, lab)
    lo_k = min(ks, key=lambda k: apm[k])
    hi_k = max(ks, key=lambda k: apm[k])
    shift = (kpm[hi_k] - kpm[lo_k]) / sdk
    print("  daily mean Kp 1950-2023: sd %.2f (n = %d days)" % (sdk, len(allkp)))
    print("  quietest group (year %d): mean Ap %.2f, mean daily Kp %.2f; stormiest group (year %d): mean Ap %.2f, mean daily Kp %.2f" % (
        lo_k, apm[lo_k], kpm[lo_k], hi_k, apm[hi_k], kpm[hi_k]))
    print("  difference = %.2f sd of daily Kp -> at 0.13-0.47%% per sd: %.2f%% to %.2f%% more deaths" % (shift, 0.13 * shift, 0.47 * shift))
    print("  the preprint's grouping: mean CDR %.3f (year %d) vs %.3f (year %d), a difference of %+.1f%%; largest spread across the 12 groups %.3f to %.3f (%+.1f%%)" % (
        mm[hi_k], hi_k, mm[lo_k], lo_k, 100.0 * (mm[hi_k] / mm[lo_k] - 1.0), min(mm.values()), max(mm.values()),
        100.0 * (max(mm.values()) / min(mm.values()) - 1.0)))
    gap = 100.0 * (mm[hi_k] / mm[lo_k] - 1.0)
    print("    that difference is %.0f to %.0f times the 263-city expectation; the 12-group spread %.0f to %.0f times" % (
        gap / (0.47 * shift), gap / (0.13 * shift),
        100.0 * (max(mm.values()) / min(mm.values()) - 1.0) / (0.47 * shift), 100.0 * (max(mm.values()) / min(mm.values()) - 1.0) / (0.13 * shift)))
    rgap = 100.0 * (rgm_cdr[hi_k] - rgm_cdr[lo_k]) / mean_rate_ns
    print("  after the trend (degree-4, fitted without the famine and COVID-19 years, which are left out): year %d minus year %d = %+.3f per 1,000, %.1f%% of the average rate, %.0f to %.0f times the 263-city expectation" % (
        hi_k, lo_k, rgm_cdr[hi_k] - rgm_cdr[lo_k], rgap, rgap / (0.47 * shift), rgap / (0.13 * shift)))

    # ------------------------------------------------------------- CSV output for the page
    if outdir:
        os.makedirs(outdir, exist_ok=True)
        with open(os.path.join(outdir, "annual_series.csv"), "w", encoding="utf-8", newline="") as f:
            w = csv.writer(f)
            w.writerow(["year", "cycle", "cycle_year_paper", "ap_mean", "cdr", "asdr_std2000", "asdr_std2023", "lex", "china_cdr"])
            for y in YEARS:
                w.writerow([y, lab[y][0], lab[y][1], "%.4f" % ap[y], "%.3f" % cdr[y], "%.4f" % asdr2000[y], "%.4f" % asdr2023[y], "%.3f" % lex[y], "%.3f" % china[y]["CDR"]])
        with open(os.path.join(outdir, "cycle_year_groups.csv"), "w", encoding="utf-8", newline="") as f:
            w = csv.writer(f)
            w.writerow(["cycle_year", "n_years", "mean_ap", "mean_cdr", "mean_asdr_std2000", "mean_lex", "mean_calendar_year", "years"])
            for k in ks:
                w.writerow([k, len(g[k]), "%.4f" % apm[k], "%.4f" % mm[k], "%.4f" % asdrm[k], "%.4f" % lexm[k], "%.2f" % yrm[k], " ".join(str(y) for y in g[k])])
        with open(os.path.join(outdir, "start_year_scan.csv"), "w", encoding="utf-8", newline="") as f:
            w = csv.writer(f)
            w.writerow(["first_year", "last_year", "r_1_12", "r_1_4", "r_5_12", "r_1_12_straight_line"])
            for row in scan_rows:
                w.writerow([row[0], row[1]] + ["%.4f" % v for v in row[2:]])
        with open(os.path.join(outdir, "null_c_r12_hist.csv"), "w", encoding="utf-8", newline="") as f:
            # counts of r(1-12) in bins 0.01 wide (lower edge inclusive), for each set of 20,000 runs; the bin
            # edges are anchored on the observed r(1-12), so the bars from the observed value up hold exactly the
            # runs that reach it
            w = csv.writer(f)
            keys = sorted(null_r12)
            anchor = r_obs["r_1_12"][0]
            w.writerow(["bin_lower", "bin_upper"] + keys)
            counts = {k: {} for k in keys}
            for k in keys:
                for v in null_r12[k]:
                    b = math.floor((v - anchor) * 100.0 + 1e-9)
                    counts[k][b] = counts[k].get(b, 0) + 1
            lo_b = min(min(c) for c in counts.values())
            hi_b = max(max(c) for c in counts.values())
            for b in range(lo_b, hi_b + 1):
                w.writerow(["%.4f" % (anchor + b / 100.0), "%.4f" % (anchor + (b + 1) / 100.0)] + [counts[k].get(b, 0) for k in keys])
        with open(os.path.join(outdir, "null_summary.csv"), "w", encoding="utf-8", newline="") as f:
            w = csv.writer(f)
            cols = ["runs", "p_r_1_12", "p_r_1_4", "p_r_5_12", "p_all_three", "median_r_1_12", "p05_r_1_12", "p95_r_1_12"]
            w.writerow(["version"] + cols)
            for k in sorted(null_summary):
                w.writerow([k] + [null_summary[k]["runs"]] + ["%.4f" % null_summary[k][c] for c in cols[1:]])
        print()
        print("CSV files written to", outdir)


if __name__ == "__main__":
    main()