#!/usr/bin/env python3
"""
The cycle test in "Is there a 6,000/12,000-year cycle?"
(Fact-Checked Science, https://factchecked.science/excursion-cycles/).

A popular claim says geomagnetic excursions come about every 6,000 years, within "a fuller harmonic up to a
12,000-year cycle". This script tests that against the dated excursions of the last 215,000 years, then looks for
any cycle length that fits better, then examines the one regular-looking stretch in the record.

What it computes, for each catalogue of excursion ages (see CATALOGUES below):
  dials        How far each event sits from the nearest mark of a cycle counted back from the present
               (marks every 6,000 years, and every 12,000 years for the "major" marks). The test statistic is
               the mean distance to the nearest mark; its chance level comes from 100,000 random catalogues of
               the same number of dates spread uniformly over 0-215,000 years. Also: how many events lie within
               1,000 years of a mark, against the number expected by chance. The claim doesn't say which marks of a
               12,000-year cycle are the major ones, so the 12,000-year dial is also scored with major marks at
               6, 18, 30, ... thousand years, and as the better of the two sets. The same tests are run on the
               best-dated events alone (one-sigma 1,000 years or less), against random dates over 0-45,000 years.
  real cycle   If a real 6,000-year cycle counted from the present had put every event on the mark nearest its age
               and only its dating error (normal, one sigma) moved it: how far from the marks the events would sit,
               and how often as far as observed (200,000 simulations), for the best-dated events and for all; the
               closest any length from 5,720 to 6,100 years comes (the video's own numbers); and how often the dial
               test above would catch such a cycle.
  hierarchy    Which events lie nearest a 12,000-year ("major") mark and which nearest a 6,000-year mark that
               is not a 12,000-year mark (a "half-cycle" mark).
  period scan  For every cycle length P from 3,000 to 60,000 years (steps of 20 years), how tightly the events
               bunch at one point of the cycle: the length R of the mean resultant of their phases 2*pi*t/P,
               with each event down-weighted for its dating uncertainty s by w = exp(-0.5*(2*pi*s/P)^2), the
               factor by which a normally distributed dating error shrinks its contribution. R runs from 0 (no
               bunching) to the mean weight (all events at the same point of the cycle), which is low for short
               cycles. So each length is judged against random catalogues at that length ("calibrated"): its local
               p is the share of random catalogues (same uncertainties) with an R at least as high there. The scan's
               statistic is the smallest local p; its chance level (p_global) comes from the same random catalogues,
               each judged the same way, and the "95% chance line" is, for each length, the R that random catalogues
               exceed so rarely that they cross the line somewhere in the scan only 5% of the time. Local p at 6,000
               and 12,000 years alone needs no choice of marks. For the checked ages, the script also simulates real
               cycles of 6,000, 6,800 and 12,000 years with random timing (1,000 each) to see how often the scan, and
               the test at that length alone, would catch them. The older "uncalibrated" numbers (one ceiling, the
               95th percentile of the best R anywhere in a scan) are kept for comparison with the original boards;
               no cycle shorter than about 10,000 years could reach that ceiling.
  recent run   The six calendar-dated events of the last 45,000 years: the spread of the gaps between them
               (coefficient of variation, CV = standard deviation / mean), against 100,000 random catalogues of
               six dates over 0-45,000 years; the same with dating errors (each age shifted 20,000 times within
               its uncertainty); and with a minimum gap between events of 3,000, 4,000 or 5,000 years.
  comb         The cycle length that best fits those six events (scan step 1 year), its teeth extended back to
               215,000 years ago and compared with the older excursions: how many teeth have an event within
               1,000 years, and how far the older events sit from the nearest tooth, against random dates.
  all events   The CV of the gaps between all events of the catalogue, against random catalogues.

Run it
  python3 cycle_test.py                 (writes cycle-test-numbers.json next to where you run it)
  python3 cycle_test.py <output.json>
  Needs Python 3 with NumPy (pip install numpy). Takes about eight minutes. Random draws use fixed seeds, so every run
  gives the same numbers.

Ages (thousands of years before AD 1950, with the one-sigma uncertainty used; sources as numbered on the page)
  "verified" catalogue, used on the page:
    Tianchi        7.0  +/- 0.45   Pan et al. 2022: 7.0 +/- 0.9 ka at 95%; possibly a local feature
    Gothenburg    14.45 +/- 0.25   Morner 1977: flip at 12,400-12,350 radiocarbon years; calibrated with IntCal20
    Hilina Pali   19.5  +/- 1.0    Scheidt et al. 2025: best-dated record 19.0-20.3 ka
    Rockall       26.5  +/- 0.5    Channell et al. 2016: 26.35-26.70 ka, one core
    Mono Lake     34.5  +/- 1.0    Liu et al. 2020; Ali et al. 2026; Nowaczyk et al. 2025
    Laschamp      41.10 +/- 0.35   Lascu et al. 2016
    Norwegian-Greenland Sea 64.5 +/- 1.5   Liu et al. 2020; Panovska et al. 2021
    Skalamaelifell 94.0 +/- 2.0    Jicha et al. 2011 (94 +/- 8 ka, 2 sigma); Shen et al. 2026 (94.47-90.95 ka)
    Post-Blake   103.5  +/- 2.0    Singer et al. 2014 (lavas about 105-104 ka); Shen et al. 2026 (103.14-99.78 ka)
    Blake        115.0  +/- 3.0    Osete et al. 2012 (116.5-112.0 ka); Shen et al. 2026 (119.1-108.8 ka)
    Iceland Basin 189.0 +/- 2.0    Channell et al. 2017 (mid-point 188.5 ka)
    Pringle Falls 211.0 +/- 6.5    Singer et al. 2008 (211 +/- 13 ka); other records give 218-239 ka
  The uncertainties of Gothenburg, Hilina Pali, Rockall, Mono Lake, the Norwegian-Greenland Sea event,
  Skalamaelifell, Post-Blake, Blake and Iceland Basin are our own estimates.
  Evidence classes (strong, real, disputed) follow the page; Iceland Basin counts as well established.
  "boards": the ages behind the original infographic boards (September 28, 2026), reproduced as a check.
  The fit-strength curves (one value per scanned cycle length) are kept for these two catalogues only.
  The other catalogues change one thing at a time (see CATALOGUES).

Convention: ages count back from AD 1950, as in the sources; the cycle's marks sit at 0, 6, 12, ... thousand
years before 1950. Counting from 2026 instead shifts every age by 76 years; the "verified_today2026" catalogue
shows that this changes nothing that matters.
"""
import json
import sys

import numpy as np

T_ALL = 215.0      # ka: the window holding all events (back to Pringle Falls)
T_YOUNG = 45.0     # ka: the window of the recent run
P_GRID = np.round(np.arange(3.0, 60.0 + 1e-9, 0.02), 4)   # cycle lengths scanned, ka

# name, age (ka before 1950), one-sigma (ka), evidence class: 'strong', 'real' (real but weaker or regional),
# 'disputed' (disputed, possibly local, or from a single record)
VERIFIED = [
    ("Tianchi", 7.0, 0.45, "disputed"),
    ("Gothenburg", 14.45, 0.25, "disputed"),
    ("Hilina Pali", 19.5, 1.0, "disputed"),
    ("Rockall", 26.5, 0.5, "disputed"),
    ("Mono Lake", 34.5, 1.0, "real"),
    ("Laschamp", 41.1, 0.35, "strong"),
    ("Norwegian-Greenland Sea", 64.5, 1.5, "real"),
    ("Skalamaelifell", 94.0, 2.0, "real"),
    ("Post-Blake", 103.5, 2.0, "real"),
    ("Blake", 115.0, 3.0, "strong"),
    ("Iceland Basin", 189.0, 2.0, "strong"),
    ("Pringle Falls", 211.0, 6.5, "real"),
]


def swap(cat, name, new):
    """Copy of a catalogue with one event replaced (new is a tuple) or removed (new is None)."""
    out = []
    for e in cat:
        if e[0] == name:
            if new is not None:
                out.append(new)
        else:
            out.append(e)
    return out


CATALOGUES = {
    # the ages used on the page
    "verified": VERIFIED,
    # the ages behind the original boards, with their evidence classes (Post-Blake as one event)
    "boards": [
        ("Tianchi", 7.5, 0.8, "disputed"), ("Gothenburg", 12.8, 0.3, "disputed"),
        ("Hilina Pali", 20.0, 1.0, "disputed"), ("Rockall", 26.5, 0.5, "disputed"),
        ("Mono Lake", 34.0, 1.5, "real"), ("Laschamp", 41.1, 0.35, "strong"),
        ("Norwegian-Greenland Sea", 64.5, 1.5, "real"), ("Post-Blake", 99.5, 5.5, "real"),
        ("Blake", 115.0, 3.0, "strong"), ("Iceland Basin", 189.0, 2.0, "real"),
        ("Pringle Falls", 211.0, 6.5, "real"),
    ],
    # Post-Blake and Skalamaelifell counted as one event spanning about 105-91 ka (as on the original boards)
    "verified_PB_merged": swap(swap(VERIFIED, "Skalamaelifell", None), "Post-Blake", ("Post-Blake", 98.0, 4.0, "real")),
    # Pringle Falls at the 238,800 years of one sediment record (Channell et al. 2012)
    "verified_PF239": swap(VERIFIED, "Pringle Falls", ("Pringle Falls", 238.8, 2.0, "real")),
    # Gothenburg at the 12,500-13,100 years of a later record that claims it (Chen et al. 2020)
    "verified_Goth12.8": swap(VERIFIED, "Gothenburg", ("Gothenburg", 12.8, 0.3, "disputed")),
    # every age counted back from 2026 instead of 1950
    "verified_today2026": [(n, round(a + 0.076, 3), s, c) for n, a, s, c in VERIFIED],
    # Hilina Pali at 20.5 ka, the middle of the 22-19 ka its records span (the page uses the best-dated record's 19.5)
    "verified_HP20.5": swap(VERIFIED, "Hilina Pali", ("Hilina Pali", 20.5, 1.0, "disputed")),
    # Blake at 114 ka, about the mid-point of the two stalagmite records (the page uses 115)
    "verified_Blake114": swap(VERIFIED, "Blake", ("Blake", 114.0, 3.0, "strong")),
}

YOUNG_NAMES = ["Tianchi", "Gothenburg", "Hilina Pali", "Rockall", "Mono Lake", "Laschamp"]

# extra subsets for the period scan: the original boards' "six better-established excursions"
EXTRA_SETS = {"boards": {"boards_six": ["Mono Lake", "Laschamp", "Norwegian-Greenland Sea", "Blake", "Iceland Basin",
                                        "Pringle Falls"]}}
# random catalogues for the scan's chance ceiling: 20,000 for the main catalogues, 4,000 for the variants
SCAN_NULL = {"verified": 20000, "boards": 20000}


def offsets(ages, P, phase=0.0):
    """Distance from each age to the nearest mark of a cycle of length P whose marks sit at phase + k*P."""
    r = np.mod(np.asarray(ages) - phase, P)
    return np.minimum(r, P - r)


def resultant(ages, sig, P):
    """Weighted mean resultant (complex) of the phases 2*pi*age/P; ages may be a 2-D array of catalogues."""
    th = 2 * np.pi * ages / P
    w = np.exp(-0.5 * (2 * np.pi * sig / P) ** 2)
    return (w * np.exp(1j * th)).sum(-1) / ages.shape[-1]


def scan(ages, sig, Ps=P_GRID):
    return np.array([abs(resultant(ages, sig, P)) for P in Ps])


def scan_max_null(n, sig, T, M, rng, Ps=P_GRID, chunk=2000):
    """Best R anywhere in the scan, for M random catalogues of n dates uniform on [0, T]."""
    out = np.empty(M)
    for i0 in range(0, M, chunk):
        sims = rng.uniform(0, T, (min(chunk, M - i0), n))
        best = np.zeros(sims.shape[0])
        for P in Ps:
            best = np.maximum(best, np.abs(resultant(sims, sig, P)))
        out[i0:i0 + sims.shape[0]] = best
    return out


def scan_null_matrix(n, sig, T, M, rng, Ps=P_GRID, chunk=2000):
    """Fit strength at every scanned length for M random catalogues of n dates uniform on [0, T] (float32, M x lengths)."""
    out = np.empty((M, len(Ps)), dtype=np.float32)
    for i0 in range(0, M, chunk):
        sims = rng.uniform(0, T, (min(chunk, M - i0), n))
        block = np.empty((sims.shape[0], len(Ps)), dtype=np.float32)
        for j, P in enumerate(Ps):
            block[:, j] = np.abs(resultant(sims, sig, P))
        out[i0:i0 + sims.shape[0]] = block
    return out


def calibrate(R, Rnull):
    """Judge each scanned length against random catalogues at that same length. The fit strength a length can reach
    depends on the length (loosely dated events count for little at short lengths), so one ceiling for all lengths
    can miss short cycles. For each length: local p = share of random catalogues at least as high there. The scan's
    statistic is the smallest local p; its chance level comes from the same random catalogues, each judged the same
    way. Sorts Rnull's columns in place, for later look-ups. Returns (summary, smallest local p of each random catalogue)."""
    M, nP = Rnull.shape
    minp_null = np.full(M, np.inf)
    p_obs = np.empty(nP)
    for j in range(nP):
        order = np.argsort(Rnull[:, j], kind="stable")
        sc = Rnull[order, j]
        rank_desc = np.empty(M)
        rank_desc[order] = np.arange(M, 0, -1)          # the highest value has rank 1
        np.minimum(minp_null, rank_desc / M, out=minp_null)
        p_obs[j] = (M - np.searchsorted(sc, np.float32(R[j]), side="left") + 1) / (M + 1)
        Rnull[:, j] = sc
    jb = int(p_obs.argmin())
    q05 = float(np.quantile(minp_null, 0.05))
    m = max(1, int(np.ceil(q05 * M)))
    line = Rnull[M - m, :].astype(float)                 # per length, the fit strength reached that rarely
    j6, j12 = int(np.argmin(np.abs(P_GRID - 6))), int(np.argmin(np.abs(P_GRID - 12)))
    cal = {
        "best_P_ka": round(float(P_GRID[jb]), 3),
        "best_local_p": round(float(p_obs[jb]), 4),
        "R_at_best": round(float(R[jb]), 4),
        "p_global": round(float((np.sum(minp_null <= p_obs[jb]) + 1) / (M + 1)), 4),
        "local_p_at_6": round(float(p_obs[j6]), 4),
        "local_p_at_12": round(float(p_obs[j12]), 4),
        "share_more_unusual_than_6": round(float((p_obs < p_obs[j6]).mean()), 4),
        "share_more_unusual_than_12": round(float((p_obs < p_obs[j12]).mean()), 4),
        "lengths_above_line": int(np.sum(R > line)),
        "local_p_threshold": round(q05, 6),
        "random_catalogues": M,
        "line_95": [round(float(x), 4) for x in line],
        "local_p": [round(float(x), 5) for x in p_obs],
    }
    return cal, minp_null


def scan_power(a, s, P0, Rnull_sorted, minp_null, rng, S=1000):
    """How often the calibrated scan would detect a real cycle of length P0 with unknown timing: every event moved to
    the cycle mark nearest its age (random phase), then shifted by a normal dating error of its one sigma."""
    M = Rnull_sorted.shape[0]
    ph = rng.uniform(0, P0, S)
    true = ph[:, None] + np.round((a[None, :] - ph[:, None]) / P0) * P0
    obs = true + rng.normal(0, 1, (S, len(a))) * s[None, :]
    minp = np.full(S, np.inf)
    j0 = int(np.argmin(np.abs(P_GRID - P0)))
    for j, P in enumerate(P_GRID):
        r = np.abs(resultant(obs, s, P)).astype(np.float32)
        pl = (M - np.searchsorted(np.ascontiguousarray(Rnull_sorted[:, j]), r, side="left") + 1) / (M + 1)
        np.minimum(minp, pl, out=minp)
        if j == j0:
            R0, p0 = r, pl
    pg = (np.searchsorted(np.sort(minp_null), minp, side="right") + 1) / (M + 1)
    return {"P_ka": P0, "share_detected_by_scan": round(float((pg < 0.05).mean()), 4),
            "share_detected_at_this_length_alone": round(float((p0 < 0.05).mean()), 4),
            "mean_R_at_P": round(float(R0.mean()), 4), "max_possible_R_at_P": round(float(np.exp(-0.5 * (2 * np.pi * s / P0) ** 2).mean()), 4),
            "simulations": S}


def real_cycle_test(a, s, P, rng, S=200000):
    """If a real cycle of length P, counted from the present, had placed every event on the mark nearest its age, and
    only dating errors (normal, one sigma) moved them: how far from the marks would they sit, and how often as far as observed?"""
    marks = np.round(a / P) * P
    obs = float(offsets(a, P).mean())
    d = offsets(marks[None, :] + rng.normal(0, 1, (S, len(a))) * s[None, :], P).mean(1)
    n_far = int(np.sum(d >= obs - 1e-12))
    return {"mean_offset_ka": round(obs, 4), "expected_if_real_ka": round(float(d.mean()), 4),
            "simulations": S, "simulations_this_far": n_far, "p_this_far_if_real": round((n_far + 1) / (S + 1), 7)}


def cv(ages):
    iv = np.diff(np.sort(np.asarray(ages, dtype=float)))
    return float(iv.std() / iv.mean()), iv


def cv_null(n, T, N, rng):
    sims = np.sort(rng.uniform(0, T, (N, n)), axis=1)
    iv = np.diff(sims, axis=1)
    return iv.std(1) / iv.mean(1)


def analyse(cat, rng, N=100000, M=20000, J=20000, extra_sets=None, keep_curve=False, power=False):
    names = [e[0] for e in cat]
    a = np.array([e[1] for e in cat], dtype=float)
    s = np.array([e[2] for e in cat], dtype=float)
    cls = [e[3] for e in cat]
    n = len(a)
    out = {"events": [{"name": nm, "age_ka": float(x), "sigma_ka": float(y), "class": c}
                      for nm, x, y, c in zip(names, a, s, cls)]}

    # ---- dials: marks every 6 and every 12 thousand years, counted from the present ----
    dials = {}
    sims = rng.uniform(0, T_ALL, (N, n))
    for P in (6.0, 12.0):
        o = offsets(a, P)
        null = offsets(sims, P).mean(1)
        within = int(np.sum(o <= 1.0 + 1e-9))
        dials[str(int(P))] = {
            "mean_offset_ka": round(float(o.mean()), 4),
            "random_mean_offset_ka": round(P / 4, 4),
            "p_random_as_close": round(float((null <= o.mean() + 1e-12).mean()), 4),
            "within_1kyr": within,
            "within_1kyr_names": [nm for nm, x in zip(names, o) if x <= 1.0 + 1e-9],
            "within_1kyr_expected_by_chance": round(n * 2.0 / P, 2),
            "per_event_offset_ka": {nm: round(float(x), 3) for nm, x in zip(names, o)},
            "per_event_phase_fraction": {nm: round(float(np.mod(x, P) / P), 4) for nm, x in zip(names, a)},
        }
    out["dials"] = dials

    # ---- hierarchy: nearest 6-kyr mark, and whether it is a 12-kyr ("major") mark ----
    hier = []
    for nm, x, c in zip(names, a, cls):
        k = x / 6.0
        lo, hi = np.floor(k) * 6.0, np.ceil(k) * 6.0
        if abs(x - lo) < abs(hi - x) - 1e-9:
            marks = [lo]
        elif abs(hi - x) < abs(x - lo) - 1e-9:
            marks = [hi]
        else:
            marks = [lo, hi]          # exactly halfway
        hier.append({"name": nm, "age_ka": float(x), "class": c,
                     "nearest_6kyr_marks_ka": [float(m) for m in marks],
                     "kind": ("halfway" if len(marks) == 2 else
                              ("major (12-kyr)" if round(marks[0]) % 12 == 0 else "half-cycle")),
                     "distance_ka": round(float(min(abs(x - m) for m in marks)), 3)})
    out["hierarchy"] = hier

    # ---- period scan, all events and the better-established ones ----
    sets = {"all": np.arange(n),
            "better_established": np.array([i for i, c in enumerate(cls) if c != "disputed"])}
    for key, members in (extra_sets or {}).items():
        sets[key] = np.array([names.index(m) for m in members])
    scans = {}
    for key, idx in sets.items():
        aa, ss = a[idx], s[idx]
        R = scan(aa, ss)
        i = int(R.argmax())
        Pb, Rb = float(P_GRID[i]), float(R[i])
        z = resultant(aa, ss, Pb)
        ph = float((np.angle(z) / (2 * np.pi) * Pb) % Pb)
        nullmax = scan_max_null(len(idx), ss, T_ALL, M, rng)
        local = np.abs(resultant(rng.uniform(0, T_ALL, (N, len(idx))), ss, Pb))
        R6 = float(R[np.argmin(np.abs(P_GRID - 6))])
        R12 = float(R[np.argmin(np.abs(P_GRID - 12))])
        lm = [j for j in range(1, len(P_GRID) - 1) if R[j] > R[j - 1] and R[j] >= R[j + 1]]
        lm = sorted(lm, key=lambda j: -R[j])[:6]
        scans[key] = {
            "names": [names[j] for j in idx],
            "best_P_ka": round(Pb, 3), "best_R": round(Rb, 4), "best_phase_ka": round(ph, 3),
            "R_at_6": round(R6, 4), "R_at_12": round(R12, 4),
            "uncalibrated_ceiling_95": round(float(np.quantile(nullmax, 0.95)), 4),
            "uncalibrated_p_global": round(float((nullmax >= Rb).mean()), 4),
            "p_local_best": round(float((local >= Rb).mean()), 4),
            "uncalibrated_share_better_than_6": round(float((R > R6).mean()), 4),
            "uncalibrated_share_better_than_12": round(float((R > R12).mean()), 4),
            "top_peaks": [(round(float(P_GRID[j]), 2), round(float(R[j]), 3)) for j in lm],
            "random_catalogues": M,
        }
        if keep_curve:
            scans[key]["curve"] = [round(float(r), 4) for r in R]
    out["scan"] = scans

    # ---- the recent run: the six calendar-dated events of the last 45 kyr ----
    yi = [names.index(nm) for nm in YOUNG_NAMES if nm in names]
    ya, ys = a[yi], s[yi]
    order = np.argsort(ya)
    ya, ys = ya[order], ys[order]
    c6, iv6 = cv(ya)
    null6 = cv_null(6, T_YOUNG, N, rng)
    aj = np.sort(ya[None, :] + rng.normal(0, 1, (J, 6)) * ys[None, :], axis=1)
    ivj = np.diff(aj, axis=1)
    cvj = ivj.std(1) / ivj.mean(1)
    q = np.quantile(cvj, [0.1, 0.5, 0.9])
    run = {"names": [YOUNG_NAMES[k] for k in order] if len(yi) == 6 else None,
           "ages_ka": [float(x) for x in ya],
           "gaps_ka": [round(float(g), 3) for g in iv6],
           "average_gap_ka": round(float(iv6.mean()), 4),
           "cv": round(c6, 4),
           "cv_random_median": round(float(np.median(null6)), 3),
           "p_cv": float((null6 <= c6).mean()),
           "cv_with_dating_errors_q10_50_90": [round(float(x), 3) for x in q],
           "p_cv_with_dating_errors_q10_50_90": [round(float((null6 <= x).mean()), 4) for x in q]}
    for d in (3, 4, 5):
        y = np.sort(rng.uniform(0, T_YOUNG - 5 * d, (N, 6)), axis=1) + np.arange(6) * d
        iv2 = np.diff(y, axis=1)
        run[f"p_cv_min_gap_{d}kyr"] = round(float(((iv2.std(1) / iv2.mean(1)) <= c6).mean()), 4)
    o6 = offsets(ya, 6.0).mean()
    n6 = offsets(rng.uniform(0, T_YOUNG, (N, 6)), 6.0).mean(1)
    run["grid6_mean_offset_ka"] = round(float(o6), 4)
    run["grid6_p_random_as_close"] = round(float((n6 <= o6 + 1e-12).mean()), 4)
    # Lake Mungo added, at its calibrated age (about 32.9 +/- 1.5 ka) and at its uncalibrated 30 ka
    for lab, age in (("LakeMungo32.9", 32.9), ("LakeMungo30.0", 30.0)):
        c7, iv7 = cv(np.append(ya, age))
        null7 = cv_null(7, T_YOUNG, N, rng)
        run[f"with_{lab}"] = {"gaps_ka": [round(float(g), 2) for g in iv7], "cv": round(c7, 3),
                              "p_cv": round(float((null7 <= c7).mean()), 4)}
    out["recent_run"] = run

    # ---- the comb fitted to the six, run back over the older events ----
    Pc = np.round(np.arange(3.0, 60.0 + 1e-9, 0.001), 4)
    Rc = np.array([abs(resultant(ya, ys, P)) for P in Pc])
    ic = int(Rc.argmax())
    P, = (float(Pc[ic]),)
    ph = float((np.angle(resultant(ya, ys, P)) / (2 * np.pi) * P) % P)
    older_idx = [i for i in range(n) if i not in yi]
    oa = a[older_idx]
    teeth = np.arange(ph, T_ALL, P)
    hits = [round(float(t), 2) for t in teeth if np.any(np.abs(a - t) <= 1.0 + 1e-9)]
    oo = offsets(oa, P, ph)
    simo = rng.uniform(T_YOUNG, T_ALL, (N, len(oa)))
    out["comb"] = {
        "period_ka": round(P, 4), "phase_ka": round(ph, 4), "R": round(float(Rc[ic]), 4),
        "young_offsets_ka": {nm: round(float(((x - ph + P / 2) % P) - P / 2), 3) for nm, x in zip(run["names"] or [], ya)},
        "teeth_0_215": int(len(teeth)), "teeth_hit_within_1kyr": len(hits), "teeth_hit_list_ka": hits,
        "teeth_hit_older_than_45ka": [t for t in hits if t > T_YOUNG],
        "older_offsets_ka": {names[i]: round(float(x), 3) for i, x in zip(older_idx, oo)},
        "older_mean_offset_ka": round(float(oo.mean()), 4),
        "chance_mean_offset_ka": round(P / 4, 4), "max_offset_ka": round(P / 2, 4),
        "p_older_this_close_by_chance": round(float((offsets(simo, P, ph).mean(1) <= oo.mean()).mean()), 4),
    }
    # could the run-back tell a cycle from none? Put each older event exactly on its nearest tooth, shift it by its
    # dating error (normal, one sigma), and score it against random dates as above
    osig = s[older_idx]
    true = ph + np.round((oa - ph) / P) * P
    nullo = np.sort(offsets(rng.uniform(T_YOUNG, T_ALL, (N, len(oa))), P, ph).mean(1))
    simt = offsets(true[None, :] + rng.normal(0, 1, (2000, len(oa))) * osig[None, :], P, ph).mean(1)
    pt = np.searchsorted(nullo, simt, side="right") / N
    out["comb"]["power"] = {"mean_offset_if_on_teeth_ka": round(float(simt.mean()), 4),
                            "share_of_simulations_with_p_below_0.05": round(float((pt < 0.05).mean()), 4),
                            "median_p_if_on_teeth": round(float(np.median(pt)), 4), "simulations": 2000}

    # ---- all events: spacing ----
    call, ivall = cv(a)
    nullall = cv_null(n, T_ALL, N, rng)
    out["all_events_spacing"] = {"gaps_ka": [round(float(g), 2) for g in ivall], "cv": round(call, 3),
                                 "cv_random_median": round(float(np.median(nullall)), 3),
                                 "p_cv": round(float((nullall <= call).mean()), 4)}

    # ---- added checks (computed last, so the random draws above are unchanged) ----
    # The claim doesn't say which marks of a 12,000-year cycle are the major ones: 12, 24, 36, ... thousand years ago
    # (above) or 6, 18, 30, ... Score the other set too (same random catalogues as the dials), and the better of the two.
    # Each event's distance to one set is 6,000 years minus its distance to the other, so the two are mirror images.
    oA = offsets(a, 12.0)
    oB = offsets(a, 12.0, 6.0)
    nA = offsets(sims, 12.0).mean(1)
    nB = offsets(sims, 12.0, 6.0).mean(1)
    out["dials"]["12"]["other_marks"] = {
        "marks_ka": "6, 18, 30, ...",
        "mean_offset_ka": round(float(oB.mean()), 4),
        "p_random_as_close": round(float((nB <= oB.mean() + 1e-12).mean()), 4),
        "within_1kyr": int(np.sum(oB <= 1.0 + 1e-9)),
        "within_1kyr_names": [nm for nm, x in zip(names, oB) if x <= 1.0 + 1e-9],
    }
    best = min(oA.mean(), oB.mean())
    out["dials"]["12"]["p_random_as_close_to_either_set"] = round(float((np.minimum(nA, nB) <= best + 1e-12).mean()), 4)
    # The best-dated events (one-sigma 1,000 years or less), against random dates over the span they cover
    wd = [i for i in range(n) if s[i] <= 1.0 + 1e-9]
    Twd = T_YOUNG if a[wd].max() <= T_YOUNG else T_ALL
    simw = rng.uniform(0, Twd, (N, len(wd)))
    bd = {"names": [names[i] for i in wd], "window_ka": Twd}
    for key, P, ph in (("6", 6.0, 0.0), ("12", 12.0, 0.0), ("12_other_marks", 12.0, 6.0)):
        o = offsets(a[wd], P, ph).mean()
        bd[key] = {"mean_offset_ka": round(float(o), 4), "random_mean_offset_ka": round(P / 4, 4),
                   "p_random_as_close": round(float((offsets(simw, P, ph).mean(1) <= o + 1e-12).mean()), 4)}
    out["dials"]["best_dated"] = bd
    # The period scan judged length by length (see calibrate), with the same number of random catalogues as above
    for key, idx in sets.items():
        aa, ss = a[idx], s[idx]
        R = scan(aa, ss)
        Rnull = scan_null_matrix(len(idx), ss, T_ALL, M, rng)
        cal, minp_null = calibrate(R, Rnull)
        if not keep_curve:
            del cal["line_95"], cal["local_p"]
        out["scan"][key]["calibrated"] = cal
        if key == "all" and power:
            out["scan"][key]["calibrated"]["power"] = [scan_power(aa, ss, P0, Rnull, minp_null, rng) for P0 in (6.0, 6.8, 12.0)]
        del Rnull
    # A real 6,000-year cycle counted from the present: the best-dated events and all of them
    rc = {"best_dated": real_cycle_test(a[wd], s[wd], 6.0, rng), "all": real_cycle_test(a, s, 6.0, rng)}
    # ... and the closest any cycle length from 5,720 to 6,100 years (the video's own numbers) comes, marks from the present
    Pv = np.round(np.arange(5.72, 6.10 + 1e-9, 0.001), 4)
    for key, idx in (("best_dated", wd), ("all", list(range(n)))):
        dd = np.array([offsets(a[idx], P).mean() for P in Pv])
        rc[key]["closest_any_length_5720_6100"] = {"mean_offset_ka": round(float(dd.min()), 4), "at_P_ka": float(Pv[dd.argmin()])}
    # how often the dial test (random dates at least as close, under 5%) would catch a real 6,000-year cycle from the present
    marks = np.round(a / 6.0) * 6.0
    dsim = offsets(marks[None, :] + rng.normal(0, 1, (20000, n)) * s[None, :], 6.0).mean(1)
    null6 = np.sort(offsets(rng.uniform(0, T_ALL, (N, n)), 6.0).mean(1))
    rc["all"]["dial_test_power"] = round(float((np.searchsorted(null6, dsim, side="right") / N < 0.05).mean()), 4)
    out["dials"]["real_cycle_6"] = rc
    return out


def summary(name, r):
    d6, d12 = r["dials"]["6"], r["dials"]["12"]
    sa, sb = r["scan"]["all"], r["scan"]["better_established"]
    rr, cb, ae = r["recent_run"], r["comb"], r["all_events_spacing"]
    L = [f"== {name}: " + ", ".join(f"{e['name']} {e['age_ka']:g}" for e in r["events"])]
    L.append(f"  6-kyr dial: mean distance {1000 * d6['mean_offset_ka']:,.0f} yr (random {1000 * d6['random_mean_offset_ka']:,.0f});"
             f" within 1 kyr {d6['within_1kyr']} of {len(r['events'])} ({', '.join(d6['within_1kyr_names'])}), chance {d6['within_1kyr_expected_by_chance']};"
             f" random dates at least as close {100 * d6['p_random_as_close']:.0f}%")
    L.append(f"  12-kyr dial: mean distance {1000 * d12['mean_offset_ka']:,.0f} yr (random {1000 * d12['random_mean_offset_ka']:,.0f});"
             f" within 1 kyr {d12['within_1kyr']} ({', '.join(d12['within_1kyr_names'])}), chance {d12['within_1kyr_expected_by_chance']};"
             f" random at least as close {100 * d12['p_random_as_close']:.0f}%")
    om, bd = d12["other_marks"], r["dials"]["best_dated"]
    L.append(f"  12-kyr dial, major marks at 6, 18, 30...: mean distance {1000 * om['mean_offset_ka']:,.0f} yr; within 1 kyr {om['within_1kyr']}"
             f" ({', '.join(om['within_1kyr_names'])}); random at least as close {100 * om['p_random_as_close']:.1f}%;"
             f" either set, the better one: {100 * d12['p_random_as_close_to_either_set']:.1f}%")
    L.append(f"  best-dated {len(bd['names'])} ({', '.join(bd['names'])}; random over 0-{bd['window_ka']:g} ka): 6-kyr {1000 * bd['6']['mean_offset_ka']:,.0f} yr,"
             f" {100 * bd['6']['p_random_as_close']:.0f}%; 12-kyr {1000 * bd['12']['mean_offset_ka']:,.0f} yr, {100 * bd['12']['p_random_as_close']:.0f}%;"
             f" 12-kyr at 6, 18, 30... {1000 * bd['12_other_marks']['mean_offset_ka']:,.0f} yr, {100 * bd['12_other_marks']['p_random_as_close']:.0f}%")
    rc = r["dials"]["real_cycle_6"]
    L.append("  a real 6-kyr cycle from the present: " + "; ".join(
        f"{k}: observed {1000 * v['mean_offset_ka']:,.0f} yr, expected {1000 * v['expected_if_real_ka']:,.0f}, this far in {v['simulations_this_far']} of"
        f" {v['simulations']:,} (p {v['p_this_far_if_real']:.6f}); closest any length 5.72-6.10: {1000 * v['closest_any_length_5720_6100']['mean_offset_ka']:,.0f} yr"
        f" at {v['closest_any_length_5720_6100']['at_P_ka']}" for k, v in rc.items() if isinstance(v, dict))
        + f"; dial test power {100 * rc['all']['dial_test_power']:.0f}%")
    L.append("  hierarchy: " + "; ".join(f"{h['name']} {h['kind']} {h['nearest_6kyr_marks_ka']}" for h in r["hierarchy"]))
    for lab, sc in [("all", sa), ("better established", sb)] + [(k, v) for k, v in r["scan"].items() if k not in ("all", "better_established")]:
        L.append(f"  scan ({lab}, n={len(sc['names'])}): best {1000 * sc['best_P_ka']:,.0f} yr, R {sc['best_R']:.3f};"
                 f" uncalibrated ceiling {sc['uncalibrated_ceiling_95']:.3f} (p {sc['uncalibrated_p_global']:.2f}); R(6) {sc['R_at_6']:.3f}, R(12) {sc['R_at_12']:.3f}")
        c = sc["calibrated"]
        L.append(f"    calibrated: best {1000 * c['best_P_ka']:,.0f} yr, local p {c['best_local_p']:.4f}, global p {c['p_global']:.2f};"
                 f" local p at 6 {c['local_p_at_6']:.2f}, at 12 {c['local_p_at_12']:.2f}; more unusual than 6: {100 * c['share_more_unusual_than_6']:.0f}%,"
                 f" than 12: {100 * c['share_more_unusual_than_12']:.0f}%; lengths above the line {c['lengths_above_line']}"
                 + "".join(f"; power at {1000 * w['P_ka']:,.0f} yr: scan {100 * w['share_detected_by_scan']:.0f}%, that length alone"
                           f" {100 * w['share_detected_at_this_length_alone']:.0f}%, mean R {w['mean_R_at_P']:.2f} (max {w['max_possible_R_at_P']:.2f})"
                           for w in c.get("power", [])))
    L.append(f"  recent run: gaps {rr['gaps_ka']} mean {1000 * rr['average_gap_ka']:,.0f}; CV {rr['cv']:.3f} (random {rr['cv_random_median']});"
             f" chance {100 * rr['p_cv']:.3f}%; with dating errors {100 * rr['p_cv_with_dating_errors_q10_50_90'][1]:.2f}%"
             f" (90th pct case {100 * rr['p_cv_with_dating_errors_q10_50_90'][2]:.1f}%); min gap 3/4/5 kyr"
             f" {100 * rr['p_cv_min_gap_3kyr']:.1f}/{100 * rr['p_cv_min_gap_4kyr']:.1f}/{100 * rr['p_cv_min_gap_5kyr']:.1f}%;"
             f" 6-kyr grid {100 * rr['grid6_p_random_as_close']:.0f}%; +Mungo 32.9 CV {rr['with_LakeMungo32.9']['cv']}, +Mungo 30.0 CV {rr['with_LakeMungo30.0']['cv']}")
    L.append(f"  comb: {1000 * cb['period_ka']:,.0f} yr, phase {cb['phase_ka']}; teeth hit {cb['teeth_hit_within_1kyr']} of {cb['teeth_0_215']}"
             f" {cb['teeth_hit_list_ka']}; older offsets {cb['older_offsets_ka']} mean {1000 * cb['older_mean_offset_ka']:,.0f}"
             f" vs chance {1000 * cb['chance_mean_offset_ka']:,.0f} (p {cb['p_older_this_close_by_chance']});"
             f" if on teeth: {1000 * cb['power']['mean_offset_if_on_teeth_ka']:,.0f}, p<0.05 in {100 * cb['power']['share_of_simulations_with_p_below_0.05']:.0f}%")
    L.append(f"  all events: CV {ae['cv']} vs random {ae['cv_random_median']} (p {ae['p_cv']})")
    return "\n".join(L)


def main(path):
    res = {"about": "Numbers behind https://factchecked.science/excursion-cycles/ (Fact-Checked Science). "
                    "Ages in thousands of years before AD 1950; see the docstring of cycle_test.py.",
           "scan_lengths_ka": [float(p) for p in P_GRID]}
    for k, (name, cat) in enumerate(CATALOGUES.items()):
        rng = np.random.default_rng(20261013 + k)     # one fixed seed per catalogue
        res[name] = analyse(cat, rng, M=SCAN_NULL.get(name, 4000), extra_sets=EXTRA_SETS.get(name), keep_curve=name in SCAN_NULL,
                            power=name == "verified")
        print(summary(name, res[name]), flush=True)
    with open(path, "w") as fh:
        json.dump(res, fh, indent=1)
    print(f"Wrote {path}")


if __name__ == "__main__":
    main(sys.argv[1] if len(sys.argv) > 1 else "cycle-test-numbers.json")