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

Our analysis for the factchecked.science page `ap-index-birth-rate`. Built from the
site's re-run of the same author's mortality preprint (`/ap-index-mortality/`,
rerun.py): the same statistics, the same year numbering and the same nulls,
applied to the world's birth rate with the preprint's 1-year lag.

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, CC BY 4.0)
    silso/TableCyclesMiMa.txt                               (WDC-SILSO, CC BY-NC 4.0; only the dates of minima are used)
    un-wpp2024/WPP2024_Demographic_Indicators_Medium.csv.gz (UN DESA, World Population Prospects 2024, CC BY 3.0 IGO)
    owid/crude-birth-rate.csv                               (optional: Our World in Data's copy of the same UN series,
                                                             which the preprint cites; checked against WPP)
    worldbank/SP.DYN.CBRT.IN_WLD.json                       (optional cross-check, including the World Bank's 2024 figure)
    un-wpp2024/WPP2024_Demographic_Indicators_OtherVariants.csv.gz
                                                            (optional: the UN's prediction intervals and other variants
                                                             for 2024-2027, against the preprint's projections)
    paper/belkin2026-ssrn-7661578.txt                       (optional: the text of the preprint's 9-page PDF, made
                                                             with `pdftotext -enc UTF-8 -layout`, to check its Tables 1 and 2;
                                                             also looked for in bryan-supplied-paper/)

With an output folder it also writes annual_series.csv, cycle_year_groups.csv,
start_year_scan.csv, lag_scan.csv, null_c_r12_hist.csv, null_summary.csv,
projections.csv and hindcast.csv there (the numbers behind the page's charts and tiles).
GFZ updates its file daily, so a fresh download can differ slightly in the latest
days from the copy we used (its checksum is in downloads/MANIFEST.md).

What it does
  1. Annual mean Ap (mean of GFZ daily Ap), 1949-2025 (2026 so far), and the UN's
     world crude birth rate (CBR, births per 1,000 people), total fertility rate
     (TFR, births per woman), births and population, 1950-2023 (estimates) and
     2024-2030 (medium projection).
  2. Numbers the years of each solar cycle from SILSO's minima under three
     conventions, groups the Ap years by number (with the next year's birth rate,
     the preprint's 1-year lag) and finds the arrangement that reproduces the
     preprint's R = 0.75 over 12 averages of 74 years; with the text of the
     preprint's PDF, it checks that arrangement against the paper's Tables 1 and 2
     (its numbering, birth rates, yearly Ap and averages) and recomputes its R.
  3. Reconstructs the preprint's projections for 2024-2027, as the paper's p. 8
     describes them, and tests the method on the years that have data (a hindcast
     from earlier base years).
  4. Tests the correlation: the trend in the birth rate, a straight line grouped
     the same way, the famine years, detrending and year-to-year changes, the TFR,
     other lags, the record's start and end, and nulls that keep the data's
     structure (shuffling whole solar cycles, shifting the solar timeline, and
     synthetic birth rates 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

AP_FIRST, AP_LAST = 1949, 2022          # the preprint's Ap years
LAG = 1                                 # births in year t paired with Ap in year t - 1
B_FIRST, B_LAST = AP_FIRST + LAG, AP_LAST + LAG   # 1950-2023, the UN's estimates
AP_YEARS = list(range(AP_FIRST, AP_LAST + 1))
B_YEARS = list(range(B_FIRST, B_LAST + 1))
PAPER_R = 0.75
PAPER_PROJ = {2024: 16.916, 2025: 17.620, 2026: 17.259, 2027: 17.501}
FAMINE = (1959, 1960, 1961)             # birth years depressed by China's famine
REBOUND = (1962, 1963)                  # the rebound that followed


# ----------------------------------------------------------------- statistics (as in the mortality re-run)

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, at=None):
    """Least-squares polynomial in scaled time; returns the fitted values (at `at` if given)."""
    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)
    pts = s if at is None else [(y - t0) / h for y in at]
    return [sum(coef[i] * si ** i for i in range(k)) for si in pts]


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, last = {}, {}, None
    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)
            last = (y, int(p[1]), int(p[2]))
    return vals, flags, last


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(path):
    world, china, india = {}, {}, {}
    keys = ("CBR", "TFR", "Births", "TPopulation1July", "MedianAgePop")
    with gzip.open(path, "rt", encoding="utf-8-sig", newline="") as f:
        for row in csv.DictReader(f):
            if row["CBR"] == "":
                continue
            t = int(row["Time"])
            if row["Location"] == "World" and row["LocTypeName"] == "World":
                world[t] = {k: float(row[k]) for k in keys}
            elif row["Location"] == "China" and row["LocTypeName"] == "Country/Area":
                china[t] = {k: float(row[k]) for k in ("CBR", "Births")}
            elif row["Location"] == "India" and row["LocTypeName"] == "Country/Area":
                india[t] = {k: float(row[k]) for k in ("CBR", "Births")}
    return world, china, india


def load_wpp_variants(path, years=(2024, 2025, 2026, 2027)):
    """World rows of WPP 2024's other projection variants: variant -> year -> (CBR, TFR, births)."""
    out = {}
    with gzip.open(path, "rt", encoding="utf-8-sig", newline="") as f:
        for row in csv.DictReader(f):
            if row["Location"] == "World" and row["LocTypeName"] == "World" and row["CBR"] != "":
                t = int(row["Time"])
                if t in years:
                    out.setdefault(row["Variant"], {})[t] = (float(row["CBR"]), float(row["TFR"]), float(row["Births"]))
    return out


def load_paper_table1(path):
    """The preprint's Table 1 (pp. 2-5 of its PDF), read from the UTF-8 text that `pdftotext -enc UTF-8 -layout` makes of the PDF:
    year -> (yearly Wolf number, ordinal number of the year in its cycle, yearly Ap, birth rate or None).
    Some Ap values are printed with a decimal comma; the table's column 6 (the birth rate a year later) is laid out
    apart from its rows and is not read here (it is column 5 moved up a row)."""
    import re
    rows = {}
    with open(path, encoding="utf-8") as f:
        for line in f.read().splitlines():  # splitlines() also splits at the form feeds between pages
            m = re.match(r"^(19[4-9]\d|20[0-2]\d)\s+([\d.]+)\s+(\d+)\s+([\d.,]+)\s+(\S+)", line)
            if m:
                b = m.group(5)
                rows[int(m.group(1))] = (float(m.group(2)), int(m.group(3)), float(m.group(4).replace(",", ".")),
                                         None if b in ("N.d.", "Н.д.") else float(b))
    return rows


def load_owid(path):
    out = {}
    with open(path, encoding="utf-8", newline="") as f:
        for row in csv.DictReader(f):
            if row["entity"] == "World":
                out[int(row["year"])] = float(row["birth_rate__sex_all__age_all__variant_estimates"])
    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 and the grouping

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 group_means(items, max_k=12):
    """items: (label k, Ap value, births value, ...) tuples -> (ks, members, mean Ap, mean births)."""
    g = {}
    for it in items:
        if it[0] <= max_k:
            g.setdefault(it[0], []).append(it)
    ks = sorted(g)
    ma = {k: sum(it[1] for it in g[k]) / len(g[k]) for k in ks}
    mb = {k: sum(it[2] for it in g[k]) / len(g[k]) for k in ks}
    return ks, g, ma, mb


def corr_sets(ks, ma, mb):
    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 ks if lo <= k <= hi]
        out[name] = (pearson([ma[k] for k in sel], [mb[k] for k in sel]), len(sel)) if len(sel) >= 3 else (float("nan"), len(sel))
    return out


def grouped_r(ap, series, lab, lag=LAG, ap_years=None, by="ap", max_k=12):
    """Group (Ap in year a, series in year a + lag) by the solar-cycle number of year a (by='ap')
    or of year a + lag (by='birth'); return the three correlations and the groups."""
    ap_years = AP_YEARS if ap_years is None else ap_years
    items = []
    for a in ap_years:
        b = a + lag
        if a in ap and b in series:
            k = lab[a][1] if by == "ap" else lab[b][1]
            items.append((k, ap[a], series[b], a, b))
    ks, g, ma, mb = group_means(items, max_k)
    return corr_sets(ks, ma, mb), ks, g, ma, mb


def show(title, r):
    print("  %-66s r(1-12) = %+.3f [n=%d, p=%s]   r(1-4) = %+.3f   r(5-12) = %+.3f" % (
        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_5_12"][0]))


def ap_years_for_births(birth_years, lag=LAG):
    return [b - lag for b in birth_years]


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

def main():
    ddir = sys.argv[1]
    outdir = sys.argv[2] if len(sys.argv) > 2 else None
    sys.stdout.reconfigure(encoding="utf-8", newline="\n")  # LF line endings on every platform

    apd, apflags, ap_last_day = 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, india = load_wpp(os.path.join(ddir, "un-wpp2024", "WPP2024_Demographic_Indicators_Medium.csv.gz"))
    owid_path = os.path.join(ddir, "owid", "crude-birth-rate.csv")
    owid = load_owid(owid_path) if os.path.exists(owid_path) else None
    paper, paper_note = None, "paper/belkin2026-ssrn-7661578.txt not found; it is optional"
    for sub in ("paper", "bryan-supplied-paper"):
        paper_path = os.path.join(ddir, sub, "belkin2026-ssrn-7661578.txt")
        if os.path.exists(paper_path):
            try:
                paper = load_paper_table1(paper_path)
            except UnicodeDecodeError:
                paper_note = "%s/belkin2026-ssrn-7661578.txt is not UTF-8 text; make it with pdftotext -enc UTF-8 -layout" % sub
            break
    wb_path = os.path.join(ddir, "worldbank", "SP.DYN.CBRT.IN_WLD.json")
    wb = load_worldbank(wb_path) if os.path.exists(wb_path) else None
    var_path = os.path.join(ddir, "un-wpp2024", "WPP2024_Demographic_Indicators_OtherVariants.csv.gz")
    variants = load_wpp_variants(var_path) if os.path.exists(var_path) else None

    ap = {y: sum(v) / len(v) for y, v in apd.items()}
    cbr = {y: wpp[y]["CBR"] for y in wpp}
    tfr = {y: wpp[y]["TFR"] for y in wpp}
    births = {y: wpp[y]["Births"] for y in wpp}
    popn = {y: wpp[y]["TPopulation1July"] for y in wpp}
    cbr_est = {y: cbr[y] for y in B_YEARS}
    tfr_est = {y: tfr[y] for y in B_YEARS}

    # ------------------------------------------------------------- 1. data
    print("=" * 110)
    print("1. DATA")
    print("=" * 110)
    print("Ap: GFZ daily Ap; last day in the file %04d-%02d-%02d; annual means %d-%d use every day (days a year %d-%d); "
          "definitive Kp on all days %d-%d: %s" % (
              ap_last_day + (AP_FIRST, 2025, min(len(apd[y]) for y in range(AP_FIRST, 2026)), max(len(apd[y]) for y in range(AP_FIRST, 2026)),
                             AP_FIRST, 2025, all(min(apflags[y]) >= 1 for y in range(AP_FIRST, 2026)))))
    print("  annual mean Ap 2019-2025: " + "  ".join("%d %.2f" % (y, ap[y]) for y in range(2019, 2026)))
    print("  2026 so far (%d days, %d definitive): mean Ap %.2f" % (len(apd[2026]), sum(1 for d in apflags[2026] if d >= 1), ap[2026]))
    print("SILSO minima (cycle, year, month):", [m for m in minima if m[0] >= 17])
    print("WPP 2024 world rows: %d-%d (estimates to %d, medium projection after)" % (min(wpp), max(wpp), B_LAST))
    if owid:
        d_ = [abs(owid[y] - cbr[y]) for y in owid if y in cbr]
        print("Our World in Data 'Birth rate' (World), %d-%d, %d years: max |OWID - WPP CBR| = %.4f per 1,000" % (
            min(owid), max(owid), len(owid), max(d_)))
    else:
        print("OWID check skipped (owid/crude-birth-rate.csv not found; it is optional)")
    if wb:
        wbd = [(y, wb[y] - cbr[y]) for y in sorted(wb) if y in cbr]
        print("World Bank SP.DYN.CBRT.IN (World) minus WPP 2024 CBR, %d-%d: max |diff| %.3f; |diff| < 0.1 in %d of %d years" % (
            wbd[0][0], wbd[-1][0], max(abs(d) for _, d in wbd), sum(1 for _, d in wbd if abs(d) < 0.1), len(wbd)))
        print("  largest differences (year, WB - WPP): " + ", ".join("%d %+.2f" % t for t in sorted(sorted(wbd, key=lambda t: -abs(t[1]))[:6])))
        print("  World Bank's latest years: " + ", ".join("%d %.3f" % (y, wb[y]) for y in sorted(wb) if y >= 2019))
    else:
        print("World Bank cross-check skipped (worldbank/SP.DYN.CBRT.IN_WLD.json not found; it is optional)")
    print()
    print("Annual series: births year t, Ap of year t-1, world CBR (per 1,000), TFR (per woman), births (millions), China CBR")
    print("  t     Ap(t-1)  CBR     TFR     births  ChinaCBR")
    for b in B_YEARS:
        print("  %d  %6.2f   %6.3f  %6.4f  %6.2f  %6.2f" % (b, ap[b - 1], cbr[b], tfr[b], births[b] / 1000.0, china[b]["CBR"]))
    print("  medium projection: " + "  ".join("%d %.3f (TFR %.3f, births %.1f M)" % (y, cbr[y], tfr[y], births[y] / 1000.0) for y in range(2024, 2028)))

    # ------------------------------------------------------------- 2. reproduce
    print()
    print("=" * 110)
    print("2. REPRODUCING R = 0.75 (12 cycle-year averages; Ap %d-%d, birth rate %d-%d)" % (AP_FIRST, AP_LAST, B_FIRST, B_LAST))
    print("=" * 110)
    labs = {c: cycle_labels(minima, c) for c in ("min", "after", "round")}
    print("  numbering  lag  grouped by  r(1-12) CBR  r(1-12) TFR   [lag L pairs Ap in year t-L with the birth rate in year t]")
    for conv in ("min", "after", "round"):
        for lag in (0, 1, 2):
            for by in ("ap", "birth"):
                ayrs = list(range(B_FIRST - lag, B_LAST - lag + 1))
                rc = grouped_r(ap, cbr_est, labs[conv], lag, ayrs, by)[0]["r_1_12"][0]
                rt = grouped_r(ap, tfr_est, labs[conv], lag, ayrs, by)[0]["r_1_12"][0]
                mark = "   <- 0.75" if round(rc, 2) == PAPER_R else ""
                print("  %-9s  %d    %-10s  %+.4f      %+.4f%s" % (conv, lag, "Ap year" if by == "ap" else "birth year", rc, rt, mark))
    lab = labs["after"]
    r_obs, ks, g, apm, bm = grouped_r(ap, cbr_est, lab)
    r_tfr, _, _, _, tm = grouped_r(ap, tfr_est, lab)
    print()
    print("The preprint's R is reproduced by numbering 'after' (year 1 = the year after the minimum year, as in the")
    print("mortality preprint), with each Ap year's number deciding the group and the next year's birth rate attached.")
    blocks = []
    cur = None
    for y in AP_YEARS:
        c = lab[y][0]
        if cur is None or cur[0] != c:
            cur = [c, []]
            blocks.append(cur)
        cur[1].append(y)
    print("Solar-cycle blocks of Ap years under this numbering:")
    for c, ys in blocks:
        print("  cycle %d: Ap years %d-%d (births %d-%d), years %d-%d of the cycle (%d years)" % (
            c, ys[0], ys[-1], ys[0] + 1, ys[-1] + 1, lab[ys[0]][1], lab[ys[-1]][1], len(ys)))
    print()
    print("The 12 data points behind R:")
    print("  yr  n  mean Ap  mean CBR  mean TFR  mean birth year  Ap years (births a year later)")
    yrm = {k: sum(it[4] for it in g[k]) / len(g[k]) for k in ks}
    for k in ks:
        print("  %2d %2d  %6.3f   %6.3f    %6.4f    %7.1f          %s" % (
            k, len(g[k]), apm[k], bm[k], tm[k], yrm[k], " ".join(str(it[3]) for it in g[k])))
    rv, n = r_obs["r_1_12"]
    print("  CBR: r(1-12) = %+.4f (paper %.2f), n = %d points from %d years, naive two-sided p = %s" % (
        rv, PAPER_R, n, sum(len(g[k]) for k in ks), fmt_p(p_two_sided(rv, n))))
    print("  CBR: r(1-4) = %+.4f (n=%d), r(5-12) = %+.4f (n=%d)  [not reported in the abstract]" % (
        r_obs["r_1_4"][0], r_obs["r_1_4"][1], r_obs["r_5_12"][0], r_obs["r_5_12"][1]))
    print("  TFR: r(1-12) = %+.4f, naive p = %s; r(1-4) = %+.4f; r(5-12) = %+.4f" % (
        r_tfr["r_1_12"][0], fmt_p(p_two_sided(*r_tfr["r_1_12"])), r_tfr["r_1_4"][0], r_tfr["r_5_12"][0]))
    print("  group means of CBR span %.3f to %.3f; of Ap %.2f to %.2f" % (min(bm.values()), max(bm.values()), min(apm.values()), max(apm.values())))

    # ------------------------------------------------------------- 2b. the paper's own tables (optional)
    print()
    print("2b. THE PAPER'S TABLES 1 AND 2 (optional: the text of its PDF, made with pdftotext -enc UTF-8 -layout)")
    k_base = lab[B_LAST - LAG][1]
    if paper:
        yrs_p = sorted(paper)
        print("  Table 1: %d rows, %d-%d" % (len(yrs_p), yrs_p[0], yrs_p[-1]))
        ordmis = [y for y in yrs_p if paper[y][1] != lab[y][1]]
        print("  its ordinal numbers against ours ('after'): %d mismatches in %d years" % (len(ordmis), len(yrs_p)))
        brmis = [y for y in yrs_p if paper[y][3] is not None and abs(paper[y][3] - cbr[y]) > 1e-9]
        print("  its birth rates against the UN's CBR: %d mismatches in %d years with data" % (
            len(brmis), sum(1 for y in yrs_p if paper[y][3] is not None)))
        dif = [(y, paper[y][2] - ap[y]) for y in yrs_p]
        big_y, big_d = max(dif, key=lambda t: abs(t[1]))
        print("  its yearly Ap minus our mean of GFZ's daily Ap: largest %+.3f (%d), mean of absolute differences %.3f; "
              "differs from ours rounded to 0.1 in %d of %d years" % (
                  big_d, big_y, statistics.mean(abs(d) for _, d in dif),
                  sum(1 for y in yrs_p if abs(round(ap[y] + 1e-9, 1) - paper[y][2]) > 1e-9), len(yrs_p)))
        ap_paper = {y: paper[y][2] for y in yrs_p}
        rp, ksp, gp, apmp, bmp = grouped_r(ap_paper, cbr_est, lab)
        print("  Table 2 rebuilt from its own Table 1: Ap means " + "  ".join("%d:%.2f" % (k, apmp[k]) for k in ksp))
        print("                                       CBR means " + "  ".join("%d:%.2f" % (k, bmp[k]) for k in ksp))
        print("  largest gap between its Ap means and ours: %.3f" % max(abs(apmp[k] - apm[k]) for k in ksp))
        print("  R with its own Ap: %+.4f (R^2 %.4f), naive two-sided p = %s; with ours: %+.4f, p = %s" % (
            rp["r_1_12"][0], rp["r_1_12"][0] ** 2, fmt_p(p_two_sided(*rp["r_1_12"])), r_obs["r_1_12"][0], fmt_p(p_two_sided(*r_obs["r_1_12"]))))
        b2 = {k: round(bm[k] + 1e-9, 2) for k in ks}
        print("  its worked example for 2024: 1.036 x 16.327 = %.4f; with Table 2's rounded averages %.2f/%.2f x 16.327 = %.4f; "
              "unrounded %.4f" % (1.036 * cbr[B_LAST], b2[4], b2[3], b2[4] / b2[3] * cbr[B_LAST], bm[4] / bm[3] * cbr[B_LAST]))
        print("  the other three from Table 2's rounded averages: " + "  ".join(
            "%d %.4f" % (t, cbr[B_LAST] * b2[lab[t - LAG][1]] / b2[k_base]) for t in (2025, 2026, 2027)))
        # the paper's method as its p. 8 states it: the ratios of successive averages (Table 2, column 6, which prints
        # three decimals), applied year after year from 2023's 16.327
        for nd in (3, 4):
            pv, chained = cbr[B_LAST], []
            for t in sorted(PAPER_PROJ):
                k = lab[t - LAG][1]
                pv *= round(bm[k] / bm[k - 1] + 1e-12, nd)
                chained.append("%d %.4f" % (t, pv))
            print("  successive ratios rounded to %d decimals, chained from 16.327: %s (paper %s)" % (
                nd, "  ".join(chained), ", ".join("%.3f" % PAPER_PROJ[t] for t in sorted(PAPER_PROJ))))
    else:
        print("  skipped (%s)" % paper_note)

    # ------------------------------------------------------------- 3. the projections
    print()
    print("=" * 110)
    print("3. THE PREPRINT'S PROJECTIONS FOR 2024-2027 (paper: 16.916, 17.620, 17.259, 17.501)")
    print("=" * 110)
    k_base = lab[B_LAST - LAG][1]
    print("Reconstruction: projected CBR(t) = CBR(2023) x [mean CBR of the group of Ap year t-1] / [mean CBR of the group of Ap year 2022]")
    print("  CBR(2023) = %.3f; Ap year 2022 is year %d of cycle %d (group mean %.4f)" % (cbr[B_LAST], k_base, lab[B_LAST - LAG][0], bm[k_base]))
    proj_rows = []
    for t in sorted(PAPER_PROJ):
        k = lab[t - LAG][1]
        v = cbr[B_LAST] * bm[k] / bm[k_base]
        imp = v * popn[t] / 1000.0
        wbv = wb.get(t) if wb else None
        print("  %d: Ap year %d is year %d (group mean %.4f): %.4f (paper %.3f, difference %+.4f); Ap of year %d = %.2f, not used; "
              "UN medium %.3f (%+.3f); implied births %.1f M vs UN %.1f M (%+.1f M); vs 2023's %.1f M %+.1f%%%s" % (
                  t, t - LAG, k, bm[k], v, PAPER_PROJ[t], v - PAPER_PROJ[t], t - LAG, ap[t - LAG], cbr[t], PAPER_PROJ[t] - cbr[t],
                  PAPER_PROJ[t] * popn[t] / 1e6, births[t] / 1000.0, (PAPER_PROJ[t] * popn[t] / 1e6) - births[t] / 1000.0,
                  births[B_LAST] / 1000.0, 100.0 * (PAPER_PROJ[t] * popn[t] / 1e6 / (births[B_LAST] / 1000.0) - 1.0),
                  ("; World Bank %.3f" % wbv) if wbv else ""))
        proj_rows.append((t, t - LAG, k, bm[k], v, PAPER_PROJ[t], cbr[t], PAPER_PROJ[t] * popn[t] / 1e6, births[t] / 1000.0, wbv))
    if variants:
        print("  The UN's other projections for the world (WPP 2024; PI = probabilistic prediction interval):")
        for t in sorted(PAPER_PROJ):
            vv = {k: variants[k][t] for k in variants if t in variants[k]}
            lo95, hi95 = vv["Lower 95 PI"][0], vv["Upper 95 PI"][0]
            cf, hi = vv["Constant fertility"], vv["High"]
            # the world TFR the preprint's value implies, interpolated between the constant-fertility and high variants
            tfr_imp = cf[1] + (PAPER_PROJ[t] - cf[0]) * (hi[1] - cf[1]) / (hi[0] - cf[0])
            print("    %d: 95%% interval %.3f-%.3f (80%% %.3f-%.3f); constant fertility %.3f; high variant %.3f; preprint %.3f, %+.3f above "
                  "the 95%% upper bound; implies a TFR of about %.2f (2023: %.2f)" % (
                      t, lo95, hi95, vv["Lower 80 PI"][0], vv["Upper 80 PI"][0], cf[0], hi[0], PAPER_PROJ[t], PAPER_PROJ[t] - hi95, tfr_imp, tfr[B_LAST]))
            proj_rows[sorted(PAPER_PROJ).index(t)] += (lo95, hi95, cf[0], hi[0])
    else:
        print("  (UN prediction intervals skipped: un-wpp2024/WPP2024_Demographic_Indicators_OtherVariants.csv.gz not found; it is optional)")
    for conv in ("min", "round"):
        _, _, _, _, bmc = grouped_r(ap, cbr_est, labs[conv])
        kb = labs[conv][B_LAST - LAG][1]
        print("  the same formula under numbering '%s' (R %+.3f): %s" % (
            conv, grouped_r(ap, cbr_est, labs[conv])[0]["r_1_12"][0],
            "  ".join("%d %.3f" % (t, cbr[B_LAST] * bmc[labs[conv][t - LAG][1]] / bmc[kb]) for t in sorted(PAPER_PROJ))))
    a_, b_ = statistics.linear_regression([apm[k] for k in ks], [bm[k] for k in ks])
    print("  Read literally, the line through the 12 averages (CBR = %.2f + %.3f x Ap) gives for births in 2023 (Ap 2022 = %.2f) "
          "%.1f per 1,000, against the UN's %.1f" % (b_, a_, ap[2022], b_ + a_ * ap[2022], cbr[2023]))
    # what the projection's ratios are made of: which groups carry the record's first and last years
    print("  Composition of the groups used (Ap years in cycle 18, i.e. births 1950-1955, and in cycle 25, births 2021-2023):")
    for k in sorted(set([k_base] + [lab[t - LAG][1] for t in PAPER_PROJ])):
        c18 = [it[3] for it in g[k] if lab[it[3]][0] == 18]
        c25 = [it[3] for it in g[k] if lab[it[3]][0] == 25]
        print("    group %d: %d years; cycle-18 member %s; cycle-25 member %s" % (
            k, len(g[k]), c18 or "none", c25 or "none"))
    # the same projection with the record cut to whole cycles, or with the trend removed first
    for title, ayrs in (("complete cycles only (Ap 1955-2019, births 1956-2020), base 2020", list(range(1955, 2020))),
                        ("record from Ap 1955 (births 1956-2023), base 2023", list(range(1955, AP_LAST + 1)))):
        _, ks2, g2, _, bm2 = grouped_r(ap, cbr_est, lab, LAG, ayrs)
        base_b = ayrs[-1] + LAG
        kb = lab[ayrs[-1]][1]
        vals = []
        for t in range(base_b + 1, base_b + 5):
            k = lab[t - LAG][1]
            vals.append("%d %.3f" % (t, cbr[base_b] * bm2[k] / bm2[kb]) if k in bm2 else "%d n/a" % t)
        print("  same method, %s: CBR(%d) = %.3f -> %s" % (title, base_b, cbr[base_b], "  ".join(vals)))
    tr4 = dict(zip(B_YEARS, polyfit_trend(B_YEARS, [cbr[b] for b in B_YEARS], 4)))
    det = {b: cbr[b] - tr4[b] for b in B_YEARS}
    _, _, _, _, dm = grouped_r(ap, det, lab)
    print("  group means of CBR minus its degree-4 trend: " + "  ".join("yr%d %+.3f" % (k, dm[k]) for k in ks))
    print("    the projection's steps from group %d with the trend removed: " % k_base + "  ".join(
        "%d %+.3f" % (t, dm[lab[t - LAG][1]] - dm[k_base]) for t in sorted(PAPER_PROJ)) + " per 1,000 (against %+.3f to %+.3f with it)" % (
        min(bm[lab[t - LAG][1]] - bm[k_base] for t in PAPER_PROJ), max(bm[lab[t - LAG][1]] - bm[k_base] for t in PAPER_PROJ)))

    # hindcast: the same method from earlier base years, scored against the UN's estimates
    print()
    print("Hindcast: the same method from each base year B = 1965-2019, using only Ap years up to B-1 and birth rates up to B,")
    print("projecting B+1 to B+4, scored against the UN's estimates; compared with 'no change' (CBR(B)) and with a straight")
    print("line fitted to the 10 years B-9..B")
    hind = []
    for B in range(1965, 2020):
        ayrs = list(range(AP_FIRST, B))
        _, ksB, gB, _, bmB = grouped_r(ap, cbr_est, lab, LAG, ayrs)
        k0 = lab[B - LAG][1]
        lin_y = list(range(B - 9, B + 1))
        for h in range(1, 5):
            t = B + h
            if t > B_LAST:
                continue
            k = lab[t - LAG][1]
            if k not in bmB or k0 not in bmB:
                continue
            fc = cbr[B] * bmB[k] / bmB[k0]
            lin = polyfit_trend(lin_y, [cbr[y] for y in lin_y], 1, at=[t])[0]
            hind.append((B, h, t, fc, cbr[B], lin, cbr[t]))
    for h in range(1, 5):
        rows = [x for x in hind if x[1] == h]
        mae_m = statistics.mean(abs(x[3] - x[6]) for x in rows)
        mae_p = statistics.mean(abs(x[4] - x[6]) for x in rows)
        mae_l = statistics.mean(abs(x[5] - x[6]) for x in rows)
        up_m = sum(1 for x in rows if x[3] > x[4])
        up_a = sum(1 for x in rows if x[6] > x[4])
        agree = sum(1 for x in rows if (x[3] - x[4]) * (x[6] - x[4]) > 0)
        worse = sum(1 for x in rows if abs(x[3] - x[6]) > abs(x[4] - x[6]))
        print("  horizon %d (n=%d): mean abs error %.3f (method) vs %.3f (no change) vs %.3f (10-year line); method worse than "
              "no change in %d of %d; method projects a rise in %d, the birth rate rose in %d; direction right in %d" % (
                  h, len(rows), mae_m, mae_p, mae_l, worse, len(rows), up_m, up_a, agree))
    for lo, hi in ((1965, 1989), (1990, 2019)):
        rows = [x for x in hind if lo <= x[0] <= hi]
        print("  base years %d-%d (n=%d): mean abs error %.3f (method) vs %.3f (no change) vs %.3f (10-year line)" % (
            lo, hi, len(rows), statistics.mean(abs(x[3] - x[6]) for x in rows), statistics.mean(abs(x[4] - x[6]) for x in rows),
            statistics.mean(abs(x[5] - x[6]) for x in rows)))
    print("  from the same position in the cycle as the preprint's base (base year B with Ap year B-1 = year %d of its cycle):" % k_base)
    for B in sorted(set(x[0] for x in hind)):
        if lab[B - LAG][1] != k_base:
            continue
        rows = [x for x in hind if x[0] == B]
        print("    base %d (CBR %.3f): " % (B, cbr[B]) + "  ".join("%d %.2f (actual %.2f)" % (x[2], x[3], x[6]) for x in rows))
    big = sorted(hind, key=lambda x: -(x[3] - x[6]))[:5]
    print("  largest over-projections: " + "; ".join("base %d -> %d: %.2f vs %.2f" % (x[0], x[2], x[3], x[6]) for x in big))
    rises = [x for x in hind if x[3] - x[4] >= 0.5]
    print("  projections of a rise of 0.5 per 1,000 or more over the base year: %d; the birth rate rose that much in %d of them" % (
        len(rises), sum(1 for x in rises if x[6] - x[4] >= 0.5)))
    print("  years 1951-2023 in which the UN's world CBR rose by 0.5 or more: %s" % (
        ", ".join(str(b) for b in B_YEARS[1:] if cbr[b] - cbr[b - 1] >= 0.5)))
    print("  last year with any rise: %d; years since 2000 with a rise: %s" % (
        max(b for b in B_YEARS[1:] if cbr[b] > cbr[b - 1]), ", ".join("%d (%+.3f)" % (b, cbr[b] - cbr[b - 1]) for b in B_YEARS if b > 2000 and cbr[b] > cbr[b - 1])))

    # ------------------------------------------------------------- 4. what drives the correlation
    print()
    print("=" * 110)
    print("4. WHAT DRIVES THE GROUPED CORRELATION")
    print("=" * 110)
    show("the birth year itself (negated), no link to Ap possible", grouped_r(ap, {b: -float(b) for b in B_YEARS}, lab)[0])
    show("a straight line falling 0.3 per 1,000 a year", grouped_r(ap, {b: 37.0 - 0.3 * (b - B_FIRST) for b in B_YEARS}, lab)[0])
    for deg in (1, 2, 3, 4):
        trend = dict(zip(B_YEARS, polyfit_trend(B_YEARS, [cbr[b] for b in B_YEARS], deg)))
        show("smooth CBR trend only (polynomial, degree %d)" % deg, grouped_r(ap, trend, lab)[0])
    for lam in (100.0, 400.0):
        trend = dict(zip(B_YEARS, hp_trend([cbr[b] for b in B_YEARS], lam)))
        show("smooth CBR trend only (Hodrick-Prescott, lambda %d)" % lam, grouped_r(ap, trend, lab)[0])
    show("world population (negated)", grouped_r(ap, {b: -popn[b] for b in B_YEARS}, lab)[0])
    show("world median age (negated)", grouped_r(ap, {b: -wpp[b]["MedianAgePop"] for b in B_YEARS}, lab)[0])
    _, _, _, _, trm = grouped_r(ap, tr4, lab)
    print("  group mean CBR minus group mean of the degree-4 trend (what the trend does not explain):")
    print("    " + "  ".join("yr%d %+.3f" % (k, bm[k] - trm[k]) for k in ks))
    print("  mean birth year of each group: " + "  ".join("yr%d %.1f" % (k, yrm[k]) for k in ks))
    sel = [k for k in ks if k <= 10]
    r10 = pearson([apm[k] for k in sel], [bm[k] for k in sel])
    print("  groups 1-10 only (drops groups 11 and 12, of 3 and 2 years): r = %+.3f (n=%d, p=%s)" % (r10, len(sel), fmt_p(p_two_sided(r10, len(sel)))))
    line = {b: 37.0 - 0.3 * (b - B_FIRST) for b in B_YEARS}
    print("  start and end of the record (births years; Ap a year earlier):")
    for b0, b1 in ((1950, 2023), (1951, 2023), (1952, 2023), (1953, 2023), (1954, 2023), (1955, 2023), (1956, 2023), (1960, 2023),
                   (1965, 2023), (1970, 2023), (1975, 2023), (1980, 2023), (1950, 2020), (1950, 2015), (1950, 2010), (1950, 2000),
                   (1950, 1986), (1987, 2023)):
        rr = grouped_r(ap, cbr_est, lab, LAG, ap_years_for_births(range(b0, b1 + 1)))[0]
        rl = grouped_r(ap, line, lab, LAG, ap_years_for_births(range(b0, b1 + 1)))[0]
        print("    births %d-%d: r(1-12) = %+.3f (straight line %+.3f)" % (b0, b1, rr["r_1_12"][0], rl["r_1_12"][0]))
    scan_rows = []
    for label_, spans in (("every first birth year 1950-1980, ending 2023", [(a, B_LAST) for a in range(1950, 1981)]),
                          ("every last birth year 1993-2023, starting 1950", [(B_FIRST, b) for b in range(1993, 2024)])):
        rs_ = [grouped_r(ap, cbr_est, lab, LAG, ap_years_for_births(range(a, b + 1)))[0]["r_1_12"][0] for a, b in spans]
        rl_ = [grouped_r(ap, line, lab, LAG, ap_years_for_births(range(a, b + 1)))[0]["r_1_12"][0] for a, b in spans]
        print("    %s (%d spans): r(1-12) %+.2f to %+.2f (negative in %d); a straight line %+.2f to %+.2f; the two curves correlate at %.3f" % (
            label_, len(spans), min(rs_), max(rs_), sum(1 for v in rs_ if v < 0), min(rl_), max(rl_), pearson(rs_, rl_)))
        if label_.startswith("every first"):
            scan_rows = [(a, b, v, w) for (a, b), v, w in zip(spans, rs_, rl_)]

    # ------------------------------------------------------------- 5. robustness
    print()
    print("=" * 110)
    print("5. ROBUSTNESS: SAME METHOD, OTHER CHOICES")
    print("=" * 110)
    show("as published: CBR, births 1950-2023", r_obs)
    show("total fertility rate (births per woman) instead", r_tfr)
    for title, bys in (("without the famine birth years 1959-61", [b for b in B_YEARS if b not in FAMINE]),
                       ("without 1959-63 (famine and rebound)", [b for b in B_YEARS if b not in FAMINE + REBOUND]),
                       ("without births 2021-23 (cycle 25)", [b for b in B_YEARS if b <= 2020]),
                       ("without births 1950-55 (cycle 18)", [b for b in B_YEARS if b >= 1956]),
                       ("complete cycles only (Ap 1955-2019, births 1956-2020)", list(range(1956, 2021))),
                       ("births 1970-2023", list(range(1970, 2024))),
                       ("first half, births 1950-1986", list(range(1950, 1987))),
                       ("second half, births 1987-2023", list(range(1987, 2024)))):
        show(title, grouped_r(ap, cbr_est, lab, LAG, ap_years_for_births(bys))[0])
        if title.startswith(("complete", "first half", "second half", "births 1970")):
            show("  ... a straight line over the same years", grouped_r(ap, line, lab, LAG, ap_years_for_births(bys))[0])
    vals = [cbr[b] for b in B_YEARS]
    det_rows = {}
    for deg in (2, 3, 4):
        res = {b: v - t for b, v, t in zip(B_YEARS, vals, polyfit_trend(B_YEARS, vals, deg))}
        det_rows["poly%d" % deg] = grouped_r(ap, res, lab)[0]
        show("CBR minus polynomial trend (degree %d)" % deg, det_rows["poly%d" % deg])
    for lam in (100.0, 400.0):
        res = {b: v - t for b, v, t in zip(B_YEARS, vals, hp_trend(vals, lam))}
        det_rows["hp%d" % lam] = grouped_r(ap, res, lab)[0]
        show("CBR minus Hodrick-Prescott trend (lambda %d)" % lam, det_rows["hp%d" % lam])
    # trend fitted with the famine and rebound years replaced by straight-line interpolation, those years left out
    shock = FAMINE + REBOUND
    sf = dict(cbr_est)
    for b in shock:
        sf[b] = cbr[1958] + (cbr[1964] - cbr[1958]) * (b - 1958) / 6.0
    for kind, tr in (("degree 4", polyfit_trend(B_YEARS, [sf[b] for b in B_YEARS], 4)), ("HP 100", hp_trend([sf[b] for b in B_YEARS], 100.0))):
        res = {b: cbr[b] - t for b, t in zip(B_YEARS, tr) if b not in shock}
        rr_ns, ks_ns, g_ns, ma_ns, mb_ns = grouped_r(ap, res, lab, LAG, ap_years_for_births(sorted(res)))
        show("CBR minus %s trend (1959-63 interpolated, then left out)" % kind, rr_ns)
        sel10 = [k for k in ks_ns if k <= 10]
        r10_ns = pearson([ma_ns[k] for k in sel10], [mb_ns[k] for k in sel10])
        parts = ["groups 1-10 only %+.3f (p %s)" % (r10_ns, fmt_p(p_two_sided(r10_ns, len(sel10))))]
        for lo, hi in ((1950, 1986), (1987, 2023), (1956, 2020)):
            kk = [b for b in sorted(res) if lo <= b <= hi]
            parts.append("births %d-%d %+.3f" % (lo, hi, grouped_r(ap, res, lab, LAG, ap_years_for_births(kk))[0]["r_1_12"][0]))
        for conv in ("min", "round"):
            parts.append("numbering '%s' %+.3f" % (conv, grouped_r(ap, res, labs[conv], LAG, ap_years_for_births(sorted(res)))[0]["r_1_12"][0]))
        for L in (0, 2):
            kk = [b for b in sorted(res) if (b - L) in ap]
            parts.append("lag %d %+.3f" % (L, grouped_r(ap, res, lab, L, [b - L for b in kk])[0]["r_1_12"][0]))
        print("    ... %s" % "; ".join(parts))
        print("    ... group means of the residuals: " + "  ".join("yr%d %+.2f" % (k, mb_ns[k]) for k in ks_ns))
    for name, series in (("CBR", cbr_est), ("TFR", tfr_est)):
        dser = {b: series[b] - series[b - 1] for b in B_YEARS[1:]}
        dap = {a: ap[a] - ap[a - 1] for a in AP_YEARS}
        show("%s changes against changes in Ap (like for like)" % name, grouped_r(dap, dser, lab, LAG, ap_years_for_births(B_YEARS[1:]))[0])
        keep = [b for b in B_YEARS[1:] if b not in shock and b - 1 not in shock]
        show("%s changes against changes in Ap, without changes into or out of 1959-63" % name, grouped_r(dap, dser, lab, LAG, ap_years_for_births(keep))[0])
        logr = {b: math.log(series[b] / series[b - 1]) for b in B_YEARS[1:]}
        show("%s year-on-year log change against changes in Ap" % name, grouped_r(dap, logr, lab, LAG, ap_years_for_births(B_YEARS[1:]))[0])
    print()
    print("Other lags, same numbering, Ap year decides the group, birth years 1950-2023 (lag L: Ap of year t-L):")
    lag_rows = []
    for L in (-2, -1, 0, 1, 2, 3, 4, 5):
        ayrs = [b - L for b in B_YEARS if (b - L) in ap and (b - L) <= 2025]
        rL = grouped_r(ap, cbr_est, lab, L, ayrs)[0]["r_1_12"][0]
        rLl = grouped_r(ap, line, lab, L, ayrs)[0]["r_1_12"][0]
        res4 = {b: cbr[b] - t for b, t in zip(B_YEARS, polyfit_trend(B_YEARS, vals, 4))}
        rLd = grouped_r(ap, res4, lab, L, ayrs)[0]["r_1_12"][0]
        lag_rows.append((L, rL, rLl, rLd))
        print("  lag %+d (Ap %d-%d): r(1-12) = %+.3f; straight line %+.3f; CBR minus degree-4 trend %+.3f" % (L, ayrs[0], ayrs[-1], rL, rLl, rLd))
    print()
    print("Other year numberings (CBR, lag 1, Ap year decides):")
    for conv in ("min", "round"):
        show("numbering '%s'" % conv, grouped_r(ap, cbr_est, labs[conv])[0])
        show("  ... a straight line", grouped_r(ap, line, labs[conv])[0])

    # ------------------------------------------------------------- 6. year by year
    print()
    print("=" * 110)
    print("6. YEAR BY YEAR (no grouping): Ap of year t-1 against the birth rate of year t, %d-%d" % (B_FIRST, B_LAST))
    print("=" * 110)
    xs = [ap[b - LAG] for b in B_YEARS]
    yby = {}
    for name, series in (("CBR", cbr_est), ("TFR", tfr_est)):
        ys = [series[b] for b in B_YEARS]
        rr = pearson(xs, ys)
        ne, r1x, r1y = n_eff(xs, ys)
        print("  Ap vs %s 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))))))
        for kind in ("poly2", "poly3", "poly4", "hp100", "hp400"):
            tr = polyfit_trend(B_YEARS, ys, int(kind[-1])) if kind.startswith("poly") else hp_trend(ys, float(kind[2:]))
            res = [v - t for v, t in zip(ys, tr)]
            r_ = pearson(xs, res)
            ne, _, _ = n_eff(xs, res)
            sel = [i for i, b in enumerate(B_YEARS) if b not in shock]
            r_ns = pearson([xs[i] for i in sel], [res[i] for i in sel])
            yby[(name, kind)] = (r_, ne)
            print("  Ap vs %s minus %-6s trend: r = %+.3f (naive p %s; effective n %.1f, p %s); without 1959-63 r = %+.3f (n = %d, p %s)" % (
                name, kind, r_, fmt_p(p_two_sided(r_, len(xs))), ne, fmt_p(p_two_sided(r_, max(3, int(ne)))), r_ns, len(sel), fmt_p(p_two_sided(r_ns, len(sel)))))
        dx = [ap[b - LAG] - ap[b - LAG - 1] for b in B_YEARS[1:]]
        dy = [series[b] - series[b - 1] for b in B_YEARS[1:]]
        rd = pearson(dx, dy)
        sel = [i for i, b in enumerate(B_YEARS[1:]) if b not in shock and b - 1 not in shock]
        rd2 = pearson([dx[i] for i in sel], [dy[i] for i in sel])
        print("  dAp vs d%s (first differences): r = %+.3f, n = %d, p = %s; without changes into or out of 1959-63 r = %+.3f (n = %d, p = %s)" % (
            name, rd, len(dx), fmt_p(p_two_sided(rd, len(dx))), rd2, len(sel), fmt_p(p_two_sided(rd2, len(sel)))))
    print("  Year by year at other lags, CBR minus HP(100) trend (1959-63 interpolated before fitting, then left out):")
    tr_sf = dict(zip(B_YEARS, hp_trend([sf[b] for b in B_YEARS], 100.0)))
    for L in (0, 1, 2, 3):
        sel = [b for b in B_YEARS if b not in shock and (b - L) in ap]
        r_ = pearson([ap[b - L] for b in sel], [cbr[b] - tr_sf[b] for b in sel])
        ne, _, _ = n_eff([ap[b - L] for b in sel], [cbr[b] - tr_sf[b] for b in sel])
        print("    lag %d: r = %+.3f, n = %d, naive p = %s, effective n = %.1f, p = %s" % (
            L, r_, len(sel), fmt_p(p_two_sided(r_, len(sel))), ne, fmt_p(p_two_sided(r_, max(3, int(ne))))))

    # ------------------------------------------------------------- 7. nulls
    print()
    print("=" * 110)
    print("7. NULLS THAT KEEP THE DATA'S STRUCTURE (numbering 'after', lag 1, births 1950-2023)")
    print("=" * 110)

    seq_lab = [lab[a][1] for a in AP_YEARS]
    seq_ap = [ap[a] for a in AP_YEARS]

    def stats_from_sequence(sl, sa, bvals):
        # sl, sa aligned with AP_YEARS; bvals aligned with B_YEARS (births a year later)
        ssum, scnt, msum = {}, {}, {}
        for i in range(len(AP_YEARS)):
            k = sl[i]
            if k > 12:
                continue
            ssum[k] = ssum.get(k, 0.0) + sa[i]
            msum[k] = msum.get(k, 0.0) + bvals[i]
            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)
        ge = sum(1 for r in results if r["r_1_12"] >= obs - 1e-12)
        r12 = sorted(r["r_1_12"] for r in results)
        print("  %s: %d alternatives; observed r(1-12) %+.3f" % (title, n, obs))
        print("    r(1-12) at least as high in %d (p = %.4f); median %+.3f, 5th-95th percentile %+.3f to %+.3f" % (
            ge, ge / n, statistics.median(r12), r12[int(0.05 * (n - 1))], r12[int(0.95 * (n - 1))]))
        return {"runs": n, "p_r_1_12": ge / 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 = [tr4[b] for b in B_YEARS]
    det4 = [cbr[b] - tr4[b] for b in B_YEARS]
    test_series = (
        ("CBR as published", [cbr[b] for b in B_YEARS], "cbr"),
        ("smooth trend only (degree-4 fit to CBR)", tr4_all, "trend"),
        ("straight falling line", [line[b] for b in B_YEARS], "line"),
        ("CBR minus degree-4 trend", det4, "detrended"),
        ("TFR as published", [tfr[b] for b in B_YEARS], "tfr"),
    )
    blk = [[(lab[a][1], ap[a]) for a 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))))
    null_summary = {}
    for title, series, key in test_series:
        obs = stats_from_sequence(seq_lab, seq_ap, series)["r_1_12"]
        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))
        null_summary["a_" + key] = 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, key in test_series[:3]:
        obs = stats_from_sequence(seq_lab, seq_ap, series)["r_1_12"]
        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))
        null_summary["a6_" + key] = summarize(title, obs, res_in)
    print("(b) Shift the solar timeline (year numbers with their Ap) against the birth rate by 0-73 years, wrapping round")
    for title, series, key in test_series:
        obs = stats_from_sequence(seq_lab, seq_ap, series)["r_1_12"]
        res_sh = []
        for k in range(len(AP_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))
        null_summary["b_" + key] = summarize(title, obs, res_sh)
    print("(c) Synthetic birth rates with CBR's trend and year-to-year noise but no solar signal (real Ap, real numbering)")
    rng = random.Random(20261030)
    obs = r_obs["r_1_12"][0]
    null_r12 = {}
    vals_b = [cbr[b] for b in B_YEARS]
    nY = len(B_YEARS)
    for deg in (3, 4):
        tr = polyfit_trend(B_YEARS, vals_b, deg)
        resid = [v - t for v, t in zip(vals_b, 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 in range(nY):
                if i:
                    e = phi * e + rng.gauss(0.0, innov)
                synth.append(tr[i] + e)
            sims.append(stats_from_sequence(seq_lab, seq_ap, synth))
        null_summary["c_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]
        sims = []
        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 = [tr[i] + e[i] for i in range(nY)]
            sims.append(stats_from_sequence(seq_lab, seq_ap, synth))
        null_summary["c_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: p(r(1-12) >= observed) %.3f to %.3f" % (
        min(v["p_r_1_12"] for k, v in null_summary.items() if k.startswith("c_")),
        max(v["p_r_1_12"] for k, v in null_summary.items() if k.startswith("c_"))))
    print("(d) Trend-free noise like the birth rate's departures from its trend, against the detrended correlation")
    obs_d = stats_from_sequence(seq_lab, seq_ap, det4)["r_1_12"]
    phi_d, sd_d = lag1(det4), statistics.pstdev(det4)
    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 in range(nY):
            if i:
                e = phi_d * e + rng.gauss(0.0, innov_d)
            synth.append(e)
        sims.append(stats_from_sequence(seq_lab, seq_ap, synth))
    null_summary["d_detrended_ar1"] = summarize("CBR minus degree-4 trend; null: AR(1) noise (lag-1 %.2f, sd %.3f), 20,000 runs" % (phi_d, sd_d), obs_d, sims)
    sims = []
    for _ in range(20000):
        e = []
        while len(e) < nY:
            s0 = rng.randrange(nY)
            e.extend(det4[(s0 + j) % nY] for j in range(5))
        sims.append(stats_from_sequence(seq_lab, seq_ap, e[:nY]))
    null_summary["d_detrended_blocks"] = summarize("CBR minus degree-4 trend; null: its own departures in 5-year blocks, 20,000 runs", obs_d, sims)
    # the same with the famine and rebound years 1959-63 interpolated before the trend is fitted and then left out
    yrs_ns = [b for b in B_YEARS if b not in shock]
    for kind, trv in (("degree 4", polyfit_trend(B_YEARS, [sf[b] for b in B_YEARS], 4)), ("HP 100", hp_trend([sf[b] for b in B_YEARS], 100.0))):
        res_ns = {b: cbr[b] - t for b, t in zip(B_YEARS, trv) if b not in shock}
        obs_ns = grouped_r(ap, res_ns, lab, LAG, ap_years_for_births(yrs_ns))[0]["r_1_12"][0]
        rv_ = [res_ns[b] for b in yrs_ns]
        phi_n, sd_n = lag1(rv_), statistics.pstdev(rv_)
        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, b in enumerate(yrs_ns):
                if i:
                    e = phi_n * e + rng.gauss(0.0, innov_n)
                synth[b] = e
            sims.append({"r_1_12": grouped_r(ap, synth, lab, LAG, ap_years_for_births(yrs_ns))[0]["r_1_12"][0]})
        null_summary["d_noshock_%s_ar1" % kind.replace(" ", "")] = summarize(
            "CBR minus %s trend, 1959-63 interpolated then left out (n = %d); null: AR(1) noise (lag-1 %.2f, sd %.3f), 20,000 runs" % (
                kind, len(yrs_ns), phi_n, sd_n), obs_ns, sims)
        sims = []
        nn = len(yrs_ns)
        for _ in range(20000):
            e = []
            while len(e) < nn:
                s0 = rng.randrange(nn)
                e.extend(rv_[(s0 + j) % nn] for j in range(5))
            sims.append({"r_1_12": grouped_r(ap, dict(zip(yrs_ns, e[:nn])), lab, LAG, ap_years_for_births(yrs_ns))[0]["r_1_12"][0]})
        null_summary["d_noshock_%s_blocks" % kind.replace(" ", "")] = summarize(
            "CBR minus %s trend, 1959-63 interpolated then left out; null: its own departures in 5-year blocks, 20,000 runs" % kind, obs_ns, sims)
    print("(e) Changes against changes, grouped (section 5's best-looking variant): the birth rate's own year-to-year changes,")
    print("    resampled in 5-year blocks, against the real changes in Ap and the real numbering")
    dap_all = {a: ap[a] - ap[a - 1] for a in AP_YEARS}
    for title, keep in (("all changes, births 1951-2023", B_YEARS[1:]),
                        ("without changes into or out of 1959-63", [b for b in B_YEARS[1:] if b not in shock and b - 1 not in shock])):
        dser = {b: cbr[b] - cbr[b - 1] for b in keep}
        obs_e = grouped_r(dap_all, dser, lab, LAG, ap_years_for_births(keep))[0]["r_1_12"][0]
        chg = [dser[b] for b in keep]
        nk = len(keep)
        sims = []
        for _ in range(20000):
            e = []
            while len(e) < nk:
                s0 = rng.randrange(nk)
                e.extend(chg[(s0 + j) % nk] for j in range(5))
            synth = dict(zip(keep, e[:nk]))
            sims.append({"r_1_12": grouped_r(dap_all, synth, lab, LAG, ap_years_for_births(keep))[0]["r_1_12"][0]})
        null_summary["e_changes_" + ("all" if len(keep) == len(B_YEARS) - 1 else "noshock")] = summarize(
            "changes, %s (%d changes), 20,000 runs" % (title, nk), obs_e, sims)
    keep_ns = [b for b in B_YEARS[1:] if b not in shock and b - 1 not in shock]
    dser_ns = {b: cbr[b] - cbr[b - 1] for b in keep_ns}
    dtfr_ns = {b: tfr[b] - tfr[b - 1] for b in keep_ns}
    print("  the same grouped changes (without 1959-63) under other choices:")
    for title, conv, L, ser, keep in (("TFR instead of CBR", "after", 1, dtfr_ns, keep_ns), ("numbering 'min'", "min", 1, dser_ns, keep_ns),
                                      ("numbering 'round'", "round", 1, dser_ns, keep_ns), ("lag 0", "after", 0, dser_ns, keep_ns),
                                      ("lag 2", "after", 2, dser_ns, keep_ns), ("lag 3", "after", 3, dser_ns, keep_ns),
                                      ("first half, births 1951-1986", "after", 1, dser_ns, [b for b in keep_ns if b <= 1986]),
                                      ("second half, births 1987-2023", "after", 1, dser_ns, [b for b in keep_ns if b >= 1987]),
                                      ("complete cycles, births 1956-2020", "after", 1, dser_ns, [b for b in keep_ns if 1956 <= b <= 2020])):
        sel = [b for b in keep if (b - L) in dap_all]
        rg = grouped_r(dap_all, ser, labs[conv], L, [b - L for b in sel])[0]["r_1_12"][0]
        ry = pearson([dap_all[b - L] for b in sel], [ser[b] for b in sel])
        print("    %-36s grouped r(1-12) = %+.3f; year by year r = %+.3f (n = %d, p = %s)" % (title, rg, ry, len(sel), fmt_p(p_two_sided(ry, len(sel)))))

    # ------------------------------------------------------------- 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="\n") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["birth_year", "ap_year", "cycle", "cycle_year_of_ap_year", "ap_mean", "cbr", "tfr", "births_thousands", "population_thousands", "china_cbr"])
            for b in B_YEARS:
                a = b - LAG
                w.writerow([b, a, lab[a][0], lab[a][1], "%.4f" % ap[a], "%.3f" % cbr[b], "%.4f" % tfr[b], "%.3f" % births[b], "%.3f" % popn[b], "%.3f" % china[b]["CBR"]])
        with open(os.path.join(outdir, "cycle_year_groups.csv"), "w", encoding="utf-8", newline="\n") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["cycle_year", "n_years", "mean_ap", "mean_cbr", "mean_tfr", "mean_birth_year", "ap_years"])
            for k in ks:
                w.writerow([k, len(g[k]), "%.4f" % apm[k], "%.4f" % bm[k], "%.4f" % tm[k], "%.2f" % yrm[k], " ".join(str(it[3]) for it in g[k])])
        with open(os.path.join(outdir, "start_year_scan.csv"), "w", encoding="utf-8", newline="\n") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["first_birth_year", "last_birth_year", "r_1_12", "r_1_12_straight_line"])
            for row in scan_rows:
                w.writerow([row[0], row[1], "%.4f" % row[2], "%.4f" % row[3]])
        with open(os.path.join(outdir, "lag_scan.csv"), "w", encoding="utf-8", newline="\n") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["lag_years", "r_1_12_cbr", "r_1_12_straight_line", "r_1_12_cbr_minus_deg4_trend"])
            for row in lag_rows:
                w.writerow([row[0]] + ["%.4f" % v for v in row[1:]])
        with open(os.path.join(outdir, "projections.csv"), "w", encoding="utf-8", newline="\n") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["year", "ap_year", "cycle_year_of_ap_year", "group_mean_cbr", "reconstructed", "preprint", "un_medium", "preprint_implied_births_millions",
                        "un_births_millions", "world_bank", "un_lower_95", "un_upper_95", "un_constant_fertility", "un_high_variant"])
            for row in proj_rows:
                extra = ["%.3f" % v for v in row[10:14]] if len(row) >= 14 else ["", "", "", ""]
                w.writerow([row[0], row[1], row[2], "%.4f" % row[3], "%.4f" % row[4], "%.3f" % row[5], "%.3f" % row[6], "%.2f" % row[7], "%.2f" % row[8],
                            "" if row[9] is None else "%.3f" % row[9]] + extra)
        with open(os.path.join(outdir, "hindcast.csv"), "w", encoding="utf-8", newline="\n") as f:
            w = csv.writer(f, lineterminator="\n")
            w.writerow(["base_year", "horizon", "year", "method", "no_change", "ten_year_line", "un_estimate"])
            for x in hind:
                w.writerow([x[0], x[1], x[2], "%.4f" % x[3], "%.3f" % x[4], "%.4f" % x[5], "%.3f" % x[6]])
        with open(os.path.join(outdir, "null_c_r12_hist.csv"), "w", encoding="utf-8", newline="\n") as f:
            # counts of r(1-12) in bins 0.01 wide (lower edge inclusive), 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, lineterminator="\n")
            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]:
                    bb = math.floor((v - anchor) * 100.0 + 1e-9)
                    counts[k][bb] = counts[k].get(bb, 0) + 1
            lo_b = min(min(c) for c in counts.values())
            hi_b = max(max(c) for c in counts.values())
            for bb in range(lo_b, hi_b + 1):
                w.writerow(["%.4f" % (anchor + bb / 100.0), "%.4f" % (anchor + (bb + 1) / 100.0)] + [counts[k].get(bb, 0) for k in keys])
        with open(os.path.join(outdir, "null_summary.csv"), "w", encoding="utf-8", newline="\n") as f:
            w = csv.writer(f, lineterminator="\n")
            cols = ["runs", "p_r_1_12", "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()