#!/usr/bin/env python3
"""
The pattern forecast in "When is the next magnetic excursion due?"
(Fact-Checked Science, https://factchecked.science/next-excursion/, section 2).

Six dated events from the last 45,000 years that have been proposed as geomagnetic excursions are
unusually evenly spaced. This script asks what that spacing says about the next one, and how far it
can be trusted:
  - average gap      (oldest - youngest) / 5
  - comb             the period P whose teeth line the six events up best. For each P from 3,000 to
                     60,000 years it measures how tightly the events bunch at one point of the cycle (the
                     length of the mean resultant of their phases, 2*pi*age/P), down-weighting each event
                     for its dating uncertainty s: weight = exp(-0.5*(2*pi*s/P)^2). The best P wins; its
                     teeth sit at phase + k*P, where phase = angle(mean resultant)*P/(2*pi).
  - straight line    an unweighted least-squares line through (event number, age) and its 95%
                     prediction interval one step ahead (Student's t with 4 degrees of freedom)
  - due by the gap   the youngest age minus the average gap
  - run back         the comb's teeth extended to 215,000 years ago and compared with the five older
                     excursions in our catalog: a tooth counts as hit if an event lies within 1,000 years,
                     and the older events' mean distance to a tooth is compared with random dates
  - evenness         the coefficient of variation (CV) of the five gaps, against 100,000 random
                     catalogues of six dates spread over 0-45,000 years; with dating errors (20,000
                     jittered copies of the ages); and with a minimum gap of 3,000-5,000 years between events
  - 6,000-year grid  how close the six sit to marks every 6,000 years counted from today, against random dates
  - odds per century the chance of an excursion starting in any given century if they come at random, for
                     several published counts (Poisson)
  - test power       whether the run-back test could detect a real cycle at all: the five older excursions
                     are placed exactly on their nearest teeth, shifted by their dating errors 2,000 times, and
                     scored against random dates as above (a separate random seed, so nothing else changes)
It runs the analysis for the verified ages used on the page, for the ages behind the original infographic
(September 30, 2026), and for single-age variants that show how much the answer depends on each date: other
published ages for Hilina Pali (17,000 to 20,500 years in steps), Gothenburg and Tianchi, and Mono Lake.

Run it
  python3 pattern.py
  Needs Python 3 with NumPy and SciPy (pip install numpy scipy). Takes about a minute. It prints the results
  and writes pattern.json; give another file name as an argument if you like: python3 pattern.py <output JSON>
  The random draws use a fixed seed, so every run gives the same numbers.

Ages (thousands of years before AD 1950, with the 1-sigma uncertainty we use; sources numbered as on the page)
  Tianchi      7.0 +/- 0.45  Pan et al. 2022 [23]: 7.0 +/- 0.9 ka, a 95% interval
  Gothenburg  14.45 +/- 0.25 Morner 1977 [24]: flip at 12,400-12,350 radiocarbon years BP; IntCal20 [27]
                             calibration of 12,350 +/- 100 gives 14,447 cal BP (95%: 14,095-14,925)
  Hilina Pali 19.5 +/- 1.0   Scheidt et al. 2025 [20]: type lavas about 19 ka, best-dated record 19.0-20.3 ka
  Rockall     26.5 +/- 0.5   Channell et al. 2016 [19]: 26.35-26.70 ka, from a single core
  Mono Lake   34.5 +/- 1.0   Liu et al. 2020 [16], Ali et al. 2026 [17], Nowaczyk et al. 2025 [28]
  Laschamp    41.10 +/- 0.35 Lascu et al. 2016 [15]
  The uncertainties for Gothenburg (+/- 100 radiocarbon years), Hilina Pali, Rockall and Mono Lake are our own
  estimates. The five older events (Norwegian-Greenland Sea, Post-Blake, Blake, Iceland Basin, Pringle Falls)
  use our catalog's ages; Blake's and Post-Blake's were not re-checked against primary sources for this page.

Convention: ages count back from AD 1950, as in the sources, so an age of t thousand years is the year
AD 1950 - 1000*t. Results are also given with "today = 2026", the convention the infographic used.
"""
import json
import sys

import numpy as np
from scipy import stats

OLDER = [  # name, age ka, 1-sigma ka: the older excursions in our catalog, used to run the comb back in time
    ("Norwegian-Greenland Sea", 64.5, 1.5),
    ("Post-Blake", 99.5, 5.5),
    ("Blake", 115.0, 3.0),
    ("Iceland Basin", 189.0, 2.0),
    ("Pringle Falls", 211.0, 6.5),
]

SCENARIOS = {
    # --- the ages behind the original infographic (September 30, 2026) ---
    "infographic": [("Tianchi", 7.5, 0.8), ("Gothenburg", 12.8, 0.3), ("Hilina Pali", 20.0, 1.0),
                    ("Rockall", 26.5, 0.5), ("Mono Lake", 34.0, 1.5), ("Laschamp", 41.1, 0.35)],
    # --- verified ages, used on the page (see the docstring) ---
    "verified": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 19.5, 1.0),
                 ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # --- single-age variants of the verified set ---
    # Hilina Pali at 17.0 ka: a figure that is really the Tianchi lava's superseded 17.1 ka age (Singer et al. 2014)
    "verified_HP17": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 17.0, 1.0),
                      ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # Hilina Pali at the type lavas' 19.3 +/- 1.6 ka (95%)
    "verified_HP19.3": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 19.3, 0.8),
                        ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # Hilina Pali at 20.5 ka, within the range of other records
    "verified_HP20.5": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 20.5, 1.0),
                        ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # Gothenburg at 12.8 ka, the age of a different record (Paleo-Danyang Lake, Chen et al. 2020)
    "verified_Goth12.8": [("Tianchi", 7.0, 0.45), ("Gothenburg", 12.8, 0.3), ("Hilina Pali", 19.5, 1.0),
                          ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # Tianchi at 7.37 ka, the charcoal age at the base of the lava (Pan et al. 2022)
    "verified_Tianchi7.4": [("Tianchi", 7.37, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 19.5, 1.0),
                            ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # Gothenburg at 13.9 ka, if the type section's radiocarbon ages needed a marine reservoir correction
    "verified_Goth13.9": [("Tianchi", 7.0, 0.45), ("Gothenburg", 13.9, 0.3), ("Hilina Pali", 19.5, 1.0),
                          ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # Mono Lake at 34.0 +/- 1.5 ka, as in the infographic
    "verified_Mono34": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 19.5, 1.0),
                        ("Rockall", 26.5, 0.5), ("Mono Lake", 34.0, 1.5), ("Laschamp", 41.1, 0.35)],
    # Added after the final check (kept last, so the random draws above are unchanged):
    # Hilina Pali across its published range, about 17 to 21.5 ka (Scheidt et al. 2025, Table 1), in steps
    "verified_HP17.5": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 17.5, 1.0),
                        ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    "verified_HP18": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 18.0, 1.0),
                      ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    "verified_HP18.5": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 18.5, 1.0),
                        ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    "verified_HP19": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 19.0, 1.0),
                      ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    "verified_HP20": [("Tianchi", 7.0, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 20.0, 1.0),
                      ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
    # Tianchi at about 8.1 ka, the age promoted before the 2022 redating (indirect dates on linked ash layers; Pan et al. 2022, p. 1)
    "verified_Tianchi8.1": [("Tianchi", 8.1, 0.45), ("Gothenburg", 14.45, 0.25), ("Hilina Pali", 19.5, 1.0),
                            ("Rockall", 26.5, 0.5), ("Mono Lake", 34.5, 1.0), ("Laschamp", 41.1, 0.35)],
}


def Rz(ages, sig, P):
    """Weighted mean resultant of the events' phases for period P (ages and P in ka)."""
    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 comb_fit(a, s, step, weighted=True, lo=3.0, hi=60.0):
    """Scan periods from lo to hi ka; return the best period, its phase (ka) and the five highest peaks."""
    Ps = np.arange(lo, hi + 1e-9, step)
    ss = s if weighted else np.zeros_like(s)
    R = np.array([abs(Rz(a, ss, P)) for P in Ps])
    i = int(R.argmax())
    P = float(Ps[i])
    z = Rz(a, ss, P)
    ph = float((np.angle(z) / (2 * np.pi) * P) % P)
    lm = [j for j in range(1, len(Ps) - 1) if R[j] > R[j - 1] and R[j] >= R[j + 1]]
    lm = sorted(lm, key=lambda j: -R[j])[:5]
    return {'P': round(P, 4), 'phase': round(ph, 4), 'R': round(float(R[i]), 4),
            'top_peaks': [(round(float(Ps[j]), 3), round(float(R[j]), 3)) for j in lm]}


def ols_pi(a):
    """Unweighted least-squares line age = b0 + b1*index (index 0 = youngest), with the 95% PI at index -1."""
    n = len(a)
    x = np.arange(n, dtype=float)
    b1, b0 = np.polyfit(x, a, 1)
    res = a - (b0 + b1 * x)
    s = np.sqrt((res ** 2).sum() / (n - 2))
    x0 = -1.0
    se = s * np.sqrt(1 + 1 / n + (x0 - x.mean()) ** 2 / ((x - x.mean()) ** 2).sum())
    t = stats.t.ppf(0.975, n - 2)
    pred = b0 + b1 * x0
    return {'slope': round(float(b1), 4), 'intercept': round(float(b0), 4), 's': round(float(s), 4),
            'pred_next_ka': round(float(pred), 4), 'halfwidth_ka': round(float(t * se), 4),
            'pi_ka': [round(float(pred - t * se), 4), round(float(pred + t * se), 4)]}


def wls_slope(a, s):
    x = np.arange(len(a), dtype=float)
    return round(float(np.polyfit(x, a, 1, w=1 / s)[0]), 4)


def offsets(ages, phase, P):
    """Distance from each age to the nearest tooth (ka)."""
    r = np.mod(ages - phase, P)
    return np.minimum(r, P - r)


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


def ad(ka, today):
    """Calendar year of an age in ka, counting back from 'today' (1950 for the sources' convention)."""
    return round(today - 1000 * ka)


def analyse(name, ev, rng, N=100000, J=20000):
    names = [e[0] for e in ev]
    a = np.array([e[1] for e in ev])
    s = np.array([e[2] for e in ev])
    order = np.argsort(a)
    a, s, names = a[order], s[order], [names[i] for i in order]
    out = {'events': [{'name': n, 'age_ka': float(x), 'sigma_ka': float(y)} for n, x, y in zip(names, a, s)]}
    gap = (a[-1] - a[0]) / 5
    c, iv = cv(a)
    out['gaps_ka'] = [round(float(g), 3) for g in iv]
    out['average_gap_ka'] = round(float(gap), 4)
    # comb fits
    out['comb_weighted_step0.02'] = comb_fit(a, s, 0.02)
    out['comb_weighted_step0.001'] = comb_fit(a, s, 0.001)
    out['comb_unweighted_step0.001'] = comb_fit(a, s, 0.001, weighted=False)
    cf = out['comb_weighted_step0.001']
    P, ph = cf['P'], cf['phase']
    out['excursions_implied_since_last_reversal'] = round(780 / P, 1)   # a cycle of period P over the 780,000 years since the last reversal
    out['last_tooth_ka'] = round(ph, 3)
    out['next_tooth_ka'] = round(ph - P, 3)
    out['next_tooth_AD_today2026'] = ad(ph - P, 2026)
    out['next_tooth_AD_BP1950'] = ad(ph - P, 1950)
    # the tooth nearest the present, whether just past or just ahead
    near = ph - round(ph / P) * P
    out['tooth_nearest_present_ka_BP'] = round(near, 3)
    out['tooth_nearest_present_AD_BP1950'] = ad(near, 1950)
    out['per_event_offset_from_comb_ka'] = {n: round(float(((x - ph + P / 2) % P) - P / 2), 3) for n, x in zip(names, a)}
    # straight line
    ol = ols_pi(a)
    out['ols'] = ol
    out['wls_slope_1_over_sigma'] = wls_slope(a, s)
    out['pi_halfwidth_years'] = round(1000 * ol['halfwidth_ka'])
    out['pi_line_own_center_AD_today2026'] = ad(ol['pred_next_ka'], 2026)
    out['pi_line_own_AD_range_today2026'] = [ad(ol['pi_ka'][1], 2026), ad(ol['pi_ka'][0], 2026)]
    out['pi_line_own_center_AD_BP1950'] = ad(ol['pred_next_ka'], 1950)
    out['pi_line_own_AD_range_BP1950'] = [ad(ol['pi_ka'][1], 1950), ad(ol['pi_ka'][0], 1950)]
    out['pi_on_comb_tooth_AD_today2026'] = [ad(ph - P + ol['halfwidth_ka'], 2026), ad(ph - P - ol['halfwidth_ka'], 2026)]
    out['pi_on_comb_tooth_AD_BP1950'] = [ad(ph - P + ol['halfwidth_ka'], 1950), ad(ph - P - ol['halfwidth_ka'], 1950)]
    # due by the average gap
    due = a[0] - gap
    out['due_by_avg_gap_ka'] = round(float(due), 4)
    out['due_by_avg_gap_AD_today2026'] = ad(due, 2026)
    out['due_by_avg_gap_AD_BP1950'] = ad(due, 1950)
    out['range_of_gaps_window_AD_today2026'] = [ad(a[0] - iv.min(), 2026), ad(a[0] - iv.max(), 2026)]
    # the comb run back through the older record
    older = np.array([o[1] for o in OLDER])
    allev = np.concatenate([a, older])
    teeth = np.arange(ph, 215, P)
    hit = [round(float(t), 2) for t in teeth if np.any(np.abs(allev - t) <= 1.0)]
    off = offsets(older, ph, P)
    out['teeth_0_215'] = int(len(teeth))
    out['teeth_hit_within_1kyr'] = len(hit)
    out['teeth_hit_list'] = hit
    out['teeth_hit_older_than_45ka'] = [t for t in hit if t > 45]
    out['older_offsets_ka'] = {o[0]: round(float(x), 3) for o, x in zip(OLDER, off)}
    out['older_mean_offset_ka'] = round(float(off.mean()), 3)
    out['chance_mean_offset_ka'] = round(P / 4, 3)
    out['max_offset_ka'] = round(P / 2, 3)
    sims = rng.uniform(45, 215, (N, len(older)))
    out['p_older_this_close_by_chance'] = round(float((offsets(sims, ph, P).mean(1) <= off.mean()).mean()), 4)
    # sensitivity: Pringle Falls at 238 ka, Blake at 116 or 112 ka
    for lab, idx, val in [('PF238', 4, 238.0), ('Blake116', 2, 116.0), ('Blake112', 2, 112.0)]:
        alt = older.copy()
        alt[idx] = val
        al = np.concatenate([a, alt])
        out[f'teeth_hit_{lab}'] = int(sum(1 for t in teeth if np.any(np.abs(al - t) <= 1.0)))
    # evenness of the spacing
    sims = np.sort(rng.uniform(0, 45, (N, 6)), axis=1)
    ivs = np.diff(sims, axis=1)
    cvs = ivs.std(1) / ivs.mean(1)
    out['cv'] = round(c, 4)
    out['cv_random_median'] = round(float(np.median(cvs)), 3)
    out['p_cv'] = float((cvs <= c).mean())
    aj = np.sort(a[None, :] + rng.normal(0, 1, (J, 6)) * s[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])
    out['cv_jitter_q10_50_90'] = [round(float(x), 3) for x in q]
    out['p_cv_at_jitter_q10_50_90'] = [round(float((cvs <= x).mean()), 4) for x in q]
    for d in (3, 4, 5):
        y = np.sort(rng.uniform(0, 45 - 5 * d, (N, 6)), axis=1) + np.arange(6) * d
        iv2 = np.diff(y, axis=1)
        c2 = iv2.std(1) / iv2.mean(1)
        out[f'p_cv_mingap_{d}'] = round(float((c2 <= c).mean()), 4)
    # a 6,000-year grid counted from today (marks at 0, 6, 12, ... ka)
    o6 = np.minimum(np.mod(a, 6), 6 - np.mod(a, 6)).mean()
    null = np.minimum(np.mod(rng.uniform(0, 45, (N, 6)), 6), 6 - np.mod(rng.uniform(0, 45, (N, 6)), 6))  # noqa: F841 (kept so the random draws match the original run)
    sims6 = rng.uniform(0, 45, (N, 6))
    n6 = np.minimum(np.mod(sims6, 6), 6 - np.mod(sims6, 6)).mean(1)
    out['grid6_mean_offset_ka'] = round(float(o6), 3)
    out['grid6_p_random_as_close'] = round(float((n6 <= o6).mean()), 3)
    # all 11 events together
    c11, _ = cv(allev)
    s11 = np.sort(rng.uniform(0, 215, (N, 11)), axis=1)
    i11 = np.diff(s11, axis=1)
    out['cv_all11'] = round(c11, 3)
    out['cv_all11_random_median'] = round(float(np.median(i11.std(1) / i11.mean(1))), 3)
    return out


def runback_power(ph, P, rng, n_sim=2000, n_null=100000):
    """Could the run-back test detect a real cycle? Put each older excursion exactly on its nearest tooth, shift it
    by its dating error (normal, one sigma), and score the mean distance to a tooth against random dates."""
    older = np.array([o[1] for o in OLDER])
    sig = np.array([o[2] for o in OLDER])
    true = ph + np.round((older - ph) / P) * P
    null = np.sort(offsets(rng.uniform(45, 215, (n_null, len(older))), ph, P).mean(1))
    obs = offsets(true[None, :] + rng.normal(0, 1, (n_sim, len(older))) * sig[None, :], ph, P).mean(1)
    p = np.searchsorted(null, obs, side='right') / n_null          # share of random catalogues at least as close
    return {'mean_offset_if_on_teeth_ka': round(float(obs.mean()), 3), 'chance_mean_offset_ka': round(P / 4, 3),
            'share_of_simulations_with_p_below_0.05': round(float((p < 0.05).mean()), 3),
            'median_p_if_on_teeth': round(float(np.median(p)), 3), 'simulations': n_sim}


def summary(name, r):
    lines = [f"== {name}: " + ", ".join(f"{e['name']} {e['age_ka']:g}" for e in r['events']) + " ka"]
    lines.append(f"  average gap {1000 * r['average_gap_ka']:,.0f} years; comb {1000 * r['comb_weighted_step0.001']['P']:,.0f} years,"
                 f" tooth nearest the present AD {r['tooth_nearest_present_AD_BP1950']}, next future tooth AD {r['next_tooth_AD_BP1950']}")
    lo, hi = r['pi_line_own_AD_range_BP1950']
    lines.append(f"  straight line points to AD {r['pi_line_own_center_AD_BP1950']}, 95% range +/-{r['pi_halfwidth_years']:,} years"
                 f" ({-lo + 1} BC to AD {hi})" if lo <= 0 else
                 f"  straight line points to AD {r['pi_line_own_center_AD_BP1950']}, 95% range +/-{r['pi_halfwidth_years']:,} years (AD {lo} to {hi})")
    lines.append(f"  due by the average gap: AD {r['due_by_avg_gap_AD_BP1950']}")
    lines.append(f"  run back to 215 ka: {r['teeth_hit_within_1kyr']} of {r['teeth_0_215']} teeth have an event within 1,000 years;"
                 f" older events {1000 * r['older_mean_offset_ka']:,.0f} years from a tooth on average vs {1000 * r['chance_mean_offset_ka']:,.0f}"
                 f" by chance (p = {r['p_older_this_close_by_chance']})")
    lines.append(f"  spacing CV {r['cv']} (random median {r['cv_random_median']}); chance this even {100 * r['p_cv']:.3f}%,"
                 f" with dating errors {100 * r['p_cv_at_jitter_q10_50_90'][1]:.2f}% (90th percentile case {100 * r['p_cv_at_jitter_q10_50_90'][2]:.1f}%),"
                 f" with a 4,000 / 5,000-year minimum gap {100 * r['p_cv_mingap_4']:.1f}% / {100 * r['p_cv_mingap_5']:.1f}%")
    lines.append(f"  6,000-year grid: random dates at least as close {100 * r['grid6_p_random_as_close']:.1f}% of the time;"
                 f" all 11 events: CV {r['cv_all11']} vs {r['cv_all11_random_median']} for random dates")
    return "\n".join(lines)


def main(path):
    rng = np.random.default_rng(20261011)
    res = {k: analyse(k, v, rng) for k, v in SCENARIOS.items()}
    # Poisson odds of an excursion starting in any given century, for several counts
    pois = {}
    for n, T in [(7, 780e3), (8, 780e3), (12, 780e3), (13, 780e3), (17, 780e3), (11, 210e3), (11, 211e3), (6, 45e3), (2, 45e3)]:
        lam = n / T
        pois[f'{n}_in_{int(T / 1000)}kyr'] = round(float(100 * (1 - np.exp(-lam * 100))), 3)
    res['poisson_percent_per_century'] = pois
    v = res['verified']['comb_weighted_step0.001']
    res['runback_power'] = runback_power(v['phase'], v['P'], np.random.default_rng(20261012))
    with open(path, 'w') as fh:
        json.dump(res, fh, indent=1)
    for k, r in res.items():
        if k not in ('poisson_percent_per_century', 'runback_power'):
            print(summary(k, r))
    print("Odds of an excursion starting in any given century (%):", pois)
    rp = res['runback_power']
    print(f"Run-back test power: if the older excursions sat on the teeth, dating errors alone would leave them"
          f" {1000 * rp['mean_offset_if_on_teeth_ka']:,.0f} years from a tooth on average (random dates: {1000 * rp['chance_mean_offset_ka']:,.0f});"
          f" the test would reach p < 0.05 in {100 * rp['share_of_simulations_with_p_below_0.05']:.0f}% of simulations")
    print(f"Wrote {path}")


if __name__ == '__main__':
    main(sys.argv[1] if len(sys.argv) > 1 else 'pattern.json')