"""Every number of ours on factchecked.science/holocene-russia-solar/, from the authors' own tables and public data.

Arslanov et al. (2026, Radiocarbon, doi:10.1017/RDC.2026.10268) date cold and warm phases in three bogs in
north-western Russia and say they are synchronous with grand solar minima and maxima. This script asks:

  1. The age model.  How well does the OxCal run printed in the paper's supplement (Figure S1, one bog, IntCal20)
     fit its own dates, and does Table S1 report that run?
  2. The match with the Sun.  Do the cold phases line up with grand minima more often than random dates would?
     Tests A-D below, with the rules fixed before they were scored; Test E checks the three bogs against each other.
  3. The plateau the paper uses as its example of a 300-year dating error, measured in IntCal20.

Inputs, all in one folder (the first six are hosted on the page; the last two are public downloads):
  tables_s1_s3.csv                    Tables S1-S3 of the 2026 supplement: each sample's depth, 14C age, calibrated age,
                                      modelled age and dT (deviation of mean annual temperature from today's, deg C)
  figS1_oxcal.csv                     the OxCal results table of the supplement's Figure S1 (pp. 7-10): each date's
                                      calibrated and modelled 95.4% range and its agreement index A
  phases_2022_table4.csv              the cold and warm phases ("limit values", cal BP) of Table 4 of the 2022 Russian
                                      version of the study (doi:10.26907/2542-064X.2022.1.135-165), whose Tables 1-3
                                      are the 2026 Tables S1-S3
  grand_minima_2022_table5.csv        the 20 grand solar minima the authors compared with (2022 Table 5)
  grand_maxima_2022_table6.csv        the 14 grand solar maxima (2022 Table 6)
  grand_minima_usoskin2007_table1.csv all 27 grand minima since 9500 BC in Usoskin, Solanki and Kovaltsov (2007),
                                      Table 1 (arXiv:0706.0385v1)
  steinhilber2012.txt                 Steinhilber et al. (2012) solar irradiance, from NOAA:
      https://www.ncei.noaa.gov/pub/data/paleo/climate_forcing/solar_variability/steinhilber2012.txt
  intcal20.14c                        the IntCal20 calibration curve (Reimer et al. 2020, CC BY 4.0):
      https://intcal.org/curves/intcal20.14c

Rules (fixed before scoring):
  Years BC/AD become cal BP as 1950 - year (BC negative); a minimum spans its centre +/- half its duration.
  Test A: a phase "matches" if a listed minimum overlaps it once the phase is widened by tol years on each side,
    tol = 0 and 250 (the authors say their phases agree "within 200-250 years"). Null 1: all phases shifted together
    by one random offset around the window, wrapping (keeps their spacing); null 2: each phase placed at random.
    Also warm phases against grand maxima, and, as a specificity check, warm phases against grand minima.
  Test B: grand-minimum centres inside cold and inside warm phases, against the share of time each covers; null: the
    whole list of minima shifted by one random offset (keeps their clustering).
  Test C: dT at the minima. Each bog's dT is interpolated linearly between samples at their modelled ages on a
    10-year grid; the short-term anomaly is dT minus its centred running mean over 1,500 years (1,000 and 3,000 as
    variants); the statistic is the mean anomaly at the minimum centres, averaged over the three bogs; null: the list
    of minima shifted around each bog's span. p = share of shifts at least as cold.
  Test D: correlation of the dT anomalies with irradiance anomalies (same high-pass), irradiance averaged over
    +/-50 years around each sample's age, samples pooled (standardized per bog); null: irradiance shifted around the
    window; p = share of shifts with r at least as large.
  Windows: 120-8,400 cal BP with the authors' 20 minima (they list none older than 6385 BC); 120-11,600 with
    Usoskin's 27; 120-9,400 for the irradiance record.
  Added after the first run (no change to Tests A-D): Test E, the pairwise correlation of the three bogs' short-term
    anomalies (null: one series shifted around the common span), the "at most" side of Test B and the share of time
    near a minimum. 20,000 random draws per test (5,000 for correlations); seed 20261026; the draws run in a fixed
    order, so the script prints the same numbers every time.
  Added after a fresh review (key "post_review", with its own seed, 20261027, so the numbers above don't move):
    Test A at margins of 0-250 years; the six cold events of Wanner et al. (2011) against the cold phases; how far the
    two cold phases that stay above today's temperature dip below their neighbouring samples; the authors' own
    two-standard-error check of the tabled ages and where Table S1's and the printed run's ages sit in each date's
    IntCal20 range; Test B with the 2022 text's own limits for three phases; Test C with Usoskin's minima in the
    same 120-8,400 window; and the gaps between the 2022 Table 4's phases and the cosmogenic ages it sets beside them.

Usage:  python3 -I holocene_checks.py <folder> [results.json]     (Python 3 standard library only; about 2 minutes)
"""
import bisect
import csv
import json
import math
import os
import random
import statistics as st
import sys

if len(sys.argv) < 2:
    sys.exit(__doc__)
IN = sys.argv[1]
OUT_JSON = sys.argv[2] if len(sys.argv) > 2 else None
rng = random.Random(20261026)
NDRAW = 20000


def read_csv(name):
    with open(os.path.join(IN, name), encoding='utf-8') as f:
        return list(csv.DictReader(f))


def events(name):
    ev = []
    for r in read_csv(name):
        c = 1950 - float(r['centre_year_BC_AD'])
        d = float(r['duration_yr'])
        ev.append((c - d / 2, c + d / 2, c))      # (young edge, old edge, centre) in cal BP
    return ev


phases = read_csv('phases_2022_table4.csv')
cold = [(float(p['young_calBP']), float(p['old_calBP'])) for p in phases if p['kind'] == 'cold']
warm = [(float(p['young_calBP']), float(p['old_calBP'])) for p in phases if p['kind'] == 'warm']
gmin5 = events('grand_minima_2022_table5.csv')
gmax6 = events('grand_maxima_2022_table6.csv')
gmin07 = events('grand_minima_usoskin2007_table1.csv')


# ---------- intervals on a circle (the window wraps around for the shift nulls) ----------
def pieces(a, b, W0, L):
    a0 = (a - W0) % L
    ln = b - a
    if ln >= L:
        return [(W0, W0 + L)]
    if a0 + ln <= L:
        return [(W0 + a0, W0 + a0 + ln)]
    return [(W0 + a0, W0 + L), (W0, W0 + a0 + ln - L)]


def overlaps(p, q):
    return p[0] <= q[1] and q[0] <= p[1]


def match_count(phs, evs, tol, W0, L):
    n = 0
    for (a, b) in phs:
        pp = pieces(a - tol, b + tol, W0, L)
        hit = False
        for (e0, e1, _) in evs:
            for q in pieces(e0, e1, W0, L):
                if any(overlaps(x, q) for x in pp):
                    hit = True
                    break
            if hit:
                break
        n += hit
    return n


def in_window(phs, W0, W1):
    return [(max(a, W0), min(b, W1)) for (a, b) in phs if b > W0 and a < W1]


def counts(xs, n):
    return [sum(1 for x in xs if x == k) for k in range(n + 1)]


def test_A(phs, evs, W0, W1, tols=(0, 250)):
    phs = in_window(phs, W0, W1)
    evs = [e for e in evs if e[1] > W0 and e[0] < W1]
    L = W1 - W0
    res = {'n_phases': len(phs), 'n_events': len(evs), 'window': [W0, W1]}
    for tol in tols:
        obs = match_count(phs, evs, tol, W0, L)
        c1 = c2 = 0
        null1 = []; null2 = []
        for _ in range(NDRAW):
            u = rng.uniform(0, L)
            m = match_count([(a + u, b + u) for (a, b) in phs], evs, tol, W0, L)
            null1.append(m)
            c1 += m >= obs
            phs2 = []
            for (a, b) in phs:
                s = rng.uniform(W0, W1)
                phs2.append((s, s + (b - a)))
            m2 = match_count(phs2, evs, tol, W0, L)
            null2.append(m2)
            c2 += m2 >= obs
        res[f'tol{tol}'] = {'observed': obs, 'p_shift': c1 / NDRAW, 'p_independent': c2 / NDRAW,
                            'mean_matches_shift_null': st.mean(null1),
                            'shift_null_counts_0_to_n': counts(null1, len(phs)),
                            'independent_null_counts_0_to_n': counts(null2, len(phs))}
    return res


def covered(phs, W0, W1):
    iv = sorted(in_window(phs, W0, W1))
    tot = 0; cur = None
    for a, b in iv:
        if cur is None or a > cur[1]:
            if cur: tot += cur[1] - cur[0]
            cur = [a, b]
        else:
            cur[1] = max(cur[1], b)
    if cur: tot += cur[1] - cur[0]
    return tot


def count_centres(phs, centres):
    return sum(1 for c in centres if any(a <= c <= b for (a, b) in phs))


def test_B(phs, evs, W0, W1):
    L = W1 - W0
    phs_w = in_window(phs, W0, W1)
    cen = [e[2] for e in evs if W0 <= e[2] <= W1]
    obs = count_centres(phs_w, cen)
    cov = covered(phs, W0, W1) / L
    ge = 0; le = 0; null = []
    for _ in range(NDRAW):
        u = rng.uniform(0, L)
        sh = [W0 + ((c - W0 + u) % L) for c in cen]
        m = count_centres(phs_w, sh)
        null.append(m); ge += m >= obs; le += m <= obs
    return {'n_centres': len(cen), 'inside': obs, 'share_of_time_covered': round(cov, 3),
            'expected_inside': round(cov * len(cen), 2), 'p_shift_at_least_as_many': ge / NDRAW,
            'p_shift_at_most_as_many': le / NDRAW, 'mean_inside_shift_null': round(st.mean(null), 2)}


def near_share(evs, W0, W1, tol):
    return covered([(e[0] - tol, e[1] + tol) for e in evs], W0, W1) / (W1 - W0)


# ---------- the bogs' dT series ----------
rows = read_csv('tables_s1_s3.csv')
bogs = {}
for r in rows:
    if r['mod'] in ('', None) or r['dT'] in ('', None):
        continue
    bogs.setdefault(r['bog'], []).append((float(r['mod']), float(r['dT'])))
for b in bogs:
    bogs[b].sort()


def grid_series(pts, step=10.0):
    ages = [p[0] for p in pts]; vals = [p[1] for p in pts]
    g0 = math.ceil(ages[0] / step) * step; g1 = math.floor(ages[-1] / step) * step
    xs = []; ys = []
    x = g0
    while x <= g1 + 1e-9:
        i = bisect.bisect_left(ages, x)
        if i == 0: y = vals[0]
        elif i >= len(ages): y = vals[-1]
        else:
            x0, x1 = ages[i - 1], ages[i]; y0, y1 = vals[i - 1], vals[i]
            y = y0 if x1 == x0 else y0 + (y1 - y0) * (x - x0) / (x1 - x0)
        xs.append(x); ys.append(y); x += step
    return xs, ys


def running_anom(xs, ys, win):
    half = int(round(win / 2 / (xs[1] - xs[0])))
    n = len(ys); pre = [0.0]
    for y in ys: pre.append(pre[-1] + y)
    an = []
    for i in range(n):
        lo = max(0, i - half); hi = min(n, i + half + 1)
        an.append(ys[i] - (pre[hi] - pre[lo]) / (hi - lo))
    return an


def value_at(xs, an, t):
    step = xs[1] - xs[0]
    i = int(round((t - xs[0]) / step))
    if 0 <= i < len(xs): return an[i]
    return None


def hist(vals, lo, hi, nb):
    w = (hi - lo) / nb
    h = [0] * nb
    for v in vals:
        k = int((v - lo) // w)
        if 0 <= k < nb: h[k] += 1
    return {'lo': lo, 'hi': hi, 'bins': nb, 'counts': h}


def test_C(evs, W0, W1, win):
    per = {}
    series = {}
    for b, pts in bogs.items():
        xs, ys = grid_series(pts)
        an = running_anom(xs, ys, win)
        lo, hi = max(W0, xs[0]), min(W1, xs[-1])
        series[b] = (xs, an, lo, hi)
    def stat(shift):
        means = []
        for b, (xs, an, lo, hi) in series.items():
            L = hi - lo
            cen = [e[2] for e in evs if lo <= e[2] <= hi]
            vals = []
            for c in cen:
                t = lo + ((c - lo + shift) % L)
                v = value_at(xs, an, t)
                if v is not None: vals.append(v)
            if vals: means.append(st.mean(vals))
        return st.mean(means)
    obs = stat(0.0)
    for b, (xs, an, lo, hi) in series.items():
        cen = [e[2] for e in evs if lo <= e[2] <= hi]
        vals = [value_at(xs, an, c) for c in cen]
        per[b] = {'n_minima': len(cen), 'mean_anomaly_at_minima': round(st.mean([v for v in vals if v is not None]), 3),
                  'sd_anomaly_all_times': round(st.pstdev([an[i] for i in range(len(xs)) if lo <= xs[i] <= hi]), 3)}
    Lref = W1 - W0
    le = 0; null = []
    for _ in range(NDRAW):
        u = rng.uniform(0, Lref)
        m = stat(u); null.append(m); le += m <= obs
    return {'window_years': win, 'observed_mean_anomaly_C': round(obs, 3), 'p_shift_at_least_as_cold': le / NDRAW,
            'null_mean': round(st.mean(null), 3), 'null_sd': round(st.pstdev(null), 3), 'per_bog': per,
            'null_histogram': hist(null, -0.5, 0.5, 40)}


def read_tsi():
    xs = []; ys = []
    started = False
    with open(os.path.join(IN, 'steinhilber2012.txt'), encoding='latin-1') as f:
        for line in f:
            s = line.split()
            if not started:
                if len(s) >= 7 and s[0] == 'Year' and s[5] == 'TSI':
                    started = True
                continue
            if len(s) >= 7:
                try:
                    xs.append(float(s[0])); ys.append(float(s[5]))
                except ValueError:
                    pass
    return xs, ys


def test_D(W0, W1, win):
    tx, ty = read_tsi()
    gx, gy = grid_series(sorted(zip(tx, ty)))
    tan = running_anom(gx, gy, win)
    L = W1 - W0
    def tsi_at(t):
        vals = [value_at(gx, tan, t + d) for d in range(-50, 51, 10)]
        vals = [v for v in vals if v is not None]
        return st.mean(vals) if vals else None
    samples = []
    for b, pts_b in bogs.items():
        xs, ys = grid_series(pts_b)
        an = running_anom(xs, ys, win)
        sel = [(a, value_at(xs, an, min(max(round(a / 10) * 10, xs[0]), xs[-1]))) for a, _ in pts_b if W0 <= a <= W1]
        sel = [(a, v) for a, v in sel if v is not None]
        m = st.mean(v for _, v in sel); s = st.pstdev([v for _, v in sel])
        samples += [(a, (v - m) / s) for a, v in sel]
    def corr(shift):
        xv = []; yv = []
        for a, v in samples:
            t = W0 + ((a - W0 + shift) % L)
            w = tsi_at(t)
            if w is not None: xv.append(w); yv.append(v)
        mx, my = st.mean(xv), st.mean(yv)
        sxy = sum((x - mx) * (y - my) for x, y in zip(xv, yv))
        sxx = sum((x - mx) ** 2 for x in xv); syy = sum((y - my) ** 2 for y in yv)
        return sxy / math.sqrt(sxx * syy), len(xv)
    obs, n = corr(0.0)
    ge = 0; null = []
    for _ in range(NDRAW // 4):
        u = rng.uniform(0, L)
        r, _ = corr(u); null.append(r); ge += r >= obs
    return {'window_years': win, 'n_samples': n, 'r_observed': round(obs, 3), 'p_shift_r_at_least': round(ge / len(null), 4),
            'null_r_sd': round(st.pstdev(null), 3)}


def test_E(win):
    ser = {}
    for b, pts in bogs.items():
        xs, ys = grid_series(pts)
        ser[b] = dict(zip(xs, running_anom(xs, ys, win)))
    names = sorted(ser)
    out = {}
    for i in range(len(names)):
        for j in range(i + 1, len(names)):
            A, B = ser[names[i]], ser[names[j]]
            common = [x for x in sorted(set(A) & set(B)) if 120 <= x <= 11600]
            a = [A[x] for x in common]; bb = [B[x] for x in common]
            def r(xv, yv):
                mx, my = st.mean(xv), st.mean(yv)
                sxy = sum((x - mx) * (y - my) for x, y in zip(xv, yv))
                return sxy / math.sqrt(sum((x - mx) ** 2 for x in xv) * sum((y - my) ** 2 for y in yv))
            obs = r(a, bb)
            n = len(common); ge = 0; null = []
            for _ in range(NDRAW // 4):
                k = rng.randrange(1, n)
                rr = r(a, bb[k:] + bb[:k]); null.append(rr); ge += rr >= obs
            out[f'{names[i]} vs {names[j]}'] = {'span': [common[0], common[-1]], 'r': round(obs, 3),
                                                'p_shift_r_at_least': round(ge / len(null), 4),
                                                'null_r_sd': round(st.pstdev(null), 3)}
    return out


# ---------- the random tests, in a fixed order ----------
res = {}
res['A_cold_vs_minima_table5'] = test_A(cold, gmin5, 120, 8400)
res['A_cold_vs_minima_usoskin2007'] = test_A(cold, gmin07, 120, 11600)
res['A_warm_vs_maxima_table6'] = test_A(warm, gmax6, 120, 8800)
res['A_warm_vs_minima_table5_specificity'] = test_A(warm, gmin5, 120, 8400)
res['B_minima_table5_in_cold'] = test_B(cold, gmin5, 120, 8400)
res['B_minima_table5_in_warm'] = test_B(warm, gmin5, 120, 8400)
res['B_minima_usoskin2007_in_cold'] = test_B(cold, gmin07, 120, 11600)
res['B_minima_usoskin2007_in_warm'] = test_B(warm, gmin07, 120, 11600)
res['B_maxima_table6_in_warm'] = test_B(warm, gmax6, 120, 8800)
res['B_maxima_table6_in_cold'] = test_B(cold, gmax6, 120, 8800)
res['share_of_window_within_tol_of_a_minimum_table5'] = {str(t): round(near_share(gmin5, 120, 8400, t), 3) for t in (0, 100, 200, 250)}
res['share_of_window_within_tol_of_a_minimum_usoskin2007'] = {str(t): round(near_share(gmin07, 120, 11600, t), 3) for t in (0, 100, 200, 250)}
res['E_bog_synchrony_win1500'] = test_E(1500)
for win in (1500, 1000, 3000):
    res[f'C_dT_at_minima_table5_win{win}'] = test_C(gmin5, 120, 8400, win)
    res[f'C_dT_at_minima_usoskin2007_win{win}'] = test_C(gmin07, 120, 11600, win)
    res[f'D_dT_vs_TSI_win{win}'] = test_D(120, 9400, win)
for win in (1000, 3000):
    res[f'E_bog_synchrony_win{win}'] = test_E(win)


# ---------- deterministic checks: the supplement's OxCal run, the tables, the plateau ----------
ox = read_csv('figS1_oxcal.csv')
oxd = {o['sample']: o for o in ox}
A = [float(o['A']) for o in ox]
nl = [r for r in rows if r['bog'] == 'Nikolsko-Lyutinskoye']
diffs = []
for r in nl:
    o = oxd[r['sample']]
    diffs.append(float(r['mod']) - (float(o['mod_from']) + float(o['mod_to'])) / 2)
ad = [abs(x) for x in diffs]
outside = sum(1 for o in ox if o['unmod_to'] and not (float(o['unmod_to']) <= (float(o['mod_from']) + float(o['mod_to'])) / 2 <= float(o['unmod_from'])))
worst = max(nl, key=lambda r: abs(float(r['mod']) - (float(oxd[r['sample']]['mod_from']) + float(oxd[r['sample']]['mod_to'])) / 2))
res['oxcal_run'] = {
    'dates': len(ox), 'A_below_60': sum(a < 60 for a in A), 'A_zero': sum(a == 0 for a in A), 'A_60_or_more': sum(a >= 60 for a in A),
    'modelled_range_width_median': st.median(float(o['mod_from']) - float(o['mod_to']) for o in ox),
    'modelled_range_width_min_max': [min(float(o['mod_from']) - float(o['mod_to']) for o in ox), max(float(o['mod_from']) - float(o['mod_to']) for o in ox)],
    'calibrated_range_width_median': st.median(float(o['unmod_from']) - float(o['unmod_to']) for o in ox if o['unmod_to']),
    'modelled_midpoint_outside_own_calibrated_range': outside,
    'tableS1_minus_run_midpoint_abs_median': st.median(ad), 'tableS1_minus_run_midpoint_abs_max': max(ad),
    'tableS1_vs_run_over_100_years': sum(x > 100 for x in ad),
    'largest_difference_sample': worst['sample'],
}
allmod = []; allcal = []; spans_all = []; per_bog = {}
for bog in ('Nikolsko-Lyutinskoye', 'Shirinsky Mokh', 'Sambalskoye'):
    r = sorted([x for x in rows if x['bog'] == bog and x['mod'] not in ('', None)], key=lambda x: (float(x['top_cm']) + float(x['bottom_cm'])) / 2)
    me = [float(x['mod_err']) for x in r]; ce = [float(x['cal_err']) for x in r]
    spans = []
    for i in range(len(r)):
        z = (float(r[i]['top_cm']) + float(r[i]['bottom_cm'])) / 2
        nb = []
        for j in (i - 1, i + 1):
            if 0 <= j < len(r):
                zj = (float(r[j]['top_cm']) + float(r[j]['bottom_cm'])) / 2
                if zj != z:
                    nb.append(abs(float(r[j]['mod']) - float(r[i]['mod'])) / abs(zj - z))
        if nb:
            spans.append(st.mean(nb) * (float(r[i]['bottom_cm']) - float(r[i]['top_cm'])))
    per_bog[bog] = {'samples_with_modelled_age': len(r), 'modelled_error_median': st.median(me),
                    'calibrated_error_median': st.median(ce), 'slice_span_median': round(st.median(spans)), 'slice_span_max': round(max(spans))}
    allmod += me; allcal += ce; spans_all += spans
res['tables_precision'] = {'per_bog': per_bog, 'modelled_error_median': st.median(allmod), 'modelled_error_30_or_less': sum(e <= 30 for e in allmod),
                           'samples': len(allmod), 'calibrated_error_median': st.median(allcal), 'slice_span_median': round(st.median(spans_all))}
curve = []
with open(os.path.join(IN, 'intcal20.14c'), encoding='utf-8') as f:
    for line in f:
        if line.startswith('#') or not line.strip():
            continue
        cal, c14 = [float(x) for x in line.split(',')[:2]]
        curve.append((cal, c14))
res['intcal20_plateau'] = {}
for lo, hi in ((2430, 2500), (2400, 2500)):
    inside = [c for c, a in curve if lo <= a <= hi]
    res['intcal20_plateau'][f'{lo}-{hi} 14C BP'] = {'cal_BP': [min(inside), max(inside)], 'span_years': max(inside) - min(inside)}
where = []
for (e0, e1, c) in gmin5:
    kinds = sorted(set(p['kind'] for p in phases if float(p['young_calBP']) <= c <= float(p['old_calBP'])))
    names = [p['label_en'] for p in phases if float(p['young_calBP']) <= c <= float(p['old_calBP'])]
    where.append({'centre_calBP': round(c), 'duration': round(e1 - e0), 'in': '; '.join(names) or '-', 'kind': '+'.join(kinds) or 'none'})
res['minima_table5_by_phase'] = where


# ---------- added after a fresh review (2026-10-26); its own random generator, so nothing above moves ----------
rng = random.Random(20261027)   # the test functions above draw from the global rng; rebinding it here leaves the frozen results as they were
post = {}

# 1. Test A at every margin from 0 to 250 years (cold phases against the 2022 version's 20 minima)
scan = test_A(cold, gmin5, 120, 8400, tols=(0, 25, 50, 75, 100, 133, 150, 200, 250))
post['A_margin_scan_table5'] = {k: {kk: vv for kk, vv in v.items() if not kk.endswith('counts_0_to_n')} if isinstance(v, dict) else v
                               for k, v in scan.items()}

# 2. The six cold events of Wanner et al. (2011), about 8200, 6300, 4700, 2700, 1550 and 550 years BP, against the cold phases
wanner = [(c, c, c) for c in (8200, 6300, 4700, 2700, 1550, 550)]
def binom_ge(n, p, k):
    return sum(math.comb(n, j) * p ** j * (1 - p) ** (n - j) for j in range(k, n + 1))
post['wanner_events'] = {}
for W1_ in (8400, 10000):
    b = test_B(cold, wanner, 120, W1_)
    b['binomial_p_at_least_as_many'] = round(binom_ge(b['n_centres'], b['share_of_time_covered'], b['inside']), 4)
    b['in_cold'] = [c for (_, _, c) in wanner if any(a <= c <= z for (a, z) in in_window(cold, 120, W1_))]
    b['in_warm'] = [c for (_, _, c) in wanner if any(a <= c <= z for (a, z) in in_window(warm, 120, W1_))]
    post['wanner_events'][f'window_120_{W1_}'] = b

# 3. The two cold phases whose dT stays above today's: the lowest value inside each, against the mean of the nearest
#    sample on each side outside it, bog by bog
dips = {}
for (a, z) in ((2410, 2820), (7190, 7340)):
    for bog, pts in bogs.items():
        inside = [(t, v) for t, v in pts if a <= t <= z]
        if not inside:
            continue
        older = [(t, v) for t, v in pts if t > z]; younger = [(t, v) for t, v in pts if t < a]
        if not older or not younger:
            continue
        nb = [min(older)[1], max(younger)[1]]
        low = min(v for _, v in inside)
        dips[f'{z}-{a} {bog}'] = {'inside': inside, 'neighbours': [min(older), max(younger)],
                                 'lowest_inside': low, 'neighbours_mean': round(st.mean(nb), 3), 'dip': round(st.mean(nb) - low, 2)}
post['cold_phase_dips'] = dips

# 4. The authors' own check (2022 p. 138: calibrated and modelled ages agree within two standard errors), and where
#    Table S1's ages and the printed run's sit against each date's own IntCal20 95.4% range
within2 = [r for r in rows if r['mod'] not in ('', None) and abs(float(r['mod']) - float(r['cal'])) <= 2 * float(r['cal_err'])]
s1_in = run_in = n_closed = 0
for o in ox:
    if not o['unmod_to']:
        continue
    n_closed += 1
    lo_, hi_ = float(o['unmod_to']), float(o['unmod_from'])
    s1v = float(next(r for r in nl if r['sample'] == o['sample'])['mod'])
    mid = (float(o['mod_from']) + float(o['mod_to'])) / 2
    s1_in += lo_ <= s1v <= hi_
    run_in += lo_ <= mid <= hi_
post['table_ages_checks'] = {'modelled_within_2_sigma_of_calibrated': len(within2),
                             'samples_with_modelled_age': sum(1 for r in rows if r['mod'] not in ('', None)),
                             'nikolsko_within_2_sigma': sum(1 for r in within2 if r['bog'] == 'Nikolsko-Lyutinskoye'),
                             'tableS1_inside_own_intcal20_range': s1_in, 'run_inside_own_intcal20_range': run_in, 'dates_with_closed_range': n_closed}

# 5. The 2022 text's own limits for three phases (4869-4340 and 1584-1210 cal BP, pp. 143 and 151; warm 6368-4960, p. 146)
cold_text = [(4340, 4869) if p == (4533, 4869) else (1210, 1584) if p == (1210, 1619) else p for p in cold]
warm_text = [(4960, 6368) if p == (4964, 6768) else p for p in warm]
post['text_limits_minima_in_cold'] = test_B(cold_text, gmin5, 120, 8400)
post['text_limits_minima_in_warm'] = test_B(warm_text, gmin5, 120, 8400)

# 6. Usoskin and colleagues' minima in the same window as the 2022 version's list (120-8,400 cal BP)
g07w = [e for e in gmin07 if 120 <= e[2] <= 8400]
for win in (1500, 1000, 3000):
    c_ = test_C(g07w, 120, 8400, win)
    c_.pop('null_histogram', None)
    post[f'C_dT_at_usoskin2007_minima_in_120_8400_win{win}'] = c_
post['usoskin2007_minima_in_120_8400'] = len(g07w)

# 7. The 2022 version's Table 4: the cosmogenic ages it pairs with each cold phase, in years before 2000 as printed,
#    counted here from 1950 (minus 50 years); how far outside the phase each lies (0 = inside)
pairs = {(1210, 1619): [1310], (2410, 2820): [2750, 2360], (4533, 4869): [4855, 4450],
         (7190, 7340): [7460, 7300, 7195]}
gaps = {}
for (a, z), ages in pairs.items():
    for b2k in ages:
        t = b2k - 50
        gaps[f'{b2k} b2k beside {z}-{a}'] = 0 if a <= t <= z else round(min(abs(t - a), abs(t - z)))
post['table4_gaps_years'] = gaps
post['table4_coldest_phase_interval'] = {'interval_b2k': [10300, 10150], 'from_1950': [10250, 10100], 'phase': [11589, 10260],
                                         'wholly_outside': True}
res['post_review'] = post

print(json.dumps(res, indent=1, ensure_ascii=False))
if OUT_JSON:
    with open(OUT_JSON, 'w', encoding='utf-8') as f:
        json.dump(res, f, indent=1, ensure_ascii=False)