"""Injection-recovery test of the Cenozoic deep-sea benthic d13C record (CENOGRID, Westerhold et al. 2020) for industrial-style carbon pulses. A pulse of M GtC of fossil organic carbon (d13C = -25 permil relative to the exogenic pool) released in << 10 kyr shifts mean-ocean d13C by A(M) = -25 * M / (M + M0), M0 = 40,000 GtC, and relaxes with e-folding time TAU_REC (carbon residence time; ~100 kyr). Bioturbation smooths the onset over TAU_MIX. The matched-filter amplitude of that template, fitted with a local linear baseline, is computed along the whole record; its robust local scatter defines the 'noise' made by natural (mostly orbital) variability. A synthetic pulse is 'flagged' when its fitted amplitude would be a > K-sigma negative outlier. Revision (referee points 1 and 5): (i) a second, deliberately different detector, a plain step statistic with no template fit (mean of the 50 kyr after onset minus mean of the 100 kyr before), is run through the same injection-recovery so that statements such as "the only 5-sigma event" can be checked against the choice of statistic; (ii) the sample-spacing coverage of CENOGRID itself (share of time lying in gaps no longer than X) is written out, because DM-SED (s05) holds no Cenozoic carbonate d13C and cannot serve as a census of geochemical sampling. Outputs: cenogrid_resolution.csv, cenogrid_scan.csv, cenogrid_events.csv, cenogrid_recovery.csv, cenogrid_step_events.csv, cenogrid_step_recovery.csv, cenogrid_spacing_coverage.csv, cenogrid_summary.json """ import os, io, json import numpy as np, pandas as pd from common import RAW, DATA, SEED, save_json M0 = 40000.0 # GtC, ocean + atmosphere + land biosphere/soils DSRC = -25.0 # permil, fossil organic carbon relative to the exogenic pool TAU_REC = 100.0 # kyr TAU_MIX = 7.0 # kyr (10 cm mixed layer at ~1.5 cm/kyr) PRE, POST = 150.0, 350.0 # fit window before / after onset, kyr K_SIGMA = 5.0 MASSES = [200, 300, 500, 700, 1000, 1500, 2000, 3000, 5000, 10000] EPOCHS = [(0, 2.58, "Quaternary"), (2.58, 5.33, "Pliocene"), (5.33, 23.03, "Miocene"), (23.03, 33.9, "Oligocene"), (33.9, 56.0, "Eocene"), (56.0, 66.0, "Paleocene")] def amp(M, m0=M0, dsrc=DSRC): M = np.asarray(M, float) return dsrc * M / (M + m0) def mass_from_amp(a, m0=M0, dsrc=DSRC): a = np.asarray(a, float) return m0 * a / (dsrc - a) def template(x, tau_rec=TAU_REC, tau_mix=TAU_MIX): """Unit-amplitude response at elapsed time x (kyr, >0 after onset).""" x = np.asarray(x, float) g = np.zeros_like(x) p = x > 0 g[p] = (np.exp(-x[p] / tau_rec) - np.exp(-x[p] / tau_mix)) * tau_rec / (tau_rec - tau_mix) return g def load(): p = os.path.join(RAW, "cenogrid", "cenogrid_S33.tab") lines = open(p, encoding="utf-8", errors="replace").read().split("\n") k = [i for i, l in enumerate(lines) if l.startswith("*/")][0] d = pd.read_csv(io.StringIO("\n".join(lines[k + 1:])), sep="\t", low_memory=False) age = d["Tuned time [Ma]"].values.astype(float) c = [c for c in d.columns if "13C" in c and "CorrAdjusted" in c][0] o = [c for c in d.columns if "18O" in c and "CorrAdjusted" in c][0] out = pd.DataFrame({"age_ma": age, "d13c": d[c].values.astype(float), "d18o": d[o].values.astype(float), "site": d["Event"].astype(str).str.replace(r"[A-Z]$", "", regex=True)}) n_all = len(out) out = out.dropna(subset=["d13c"]).sort_values("age_ma").reset_index(drop=True) return out, n_all def scan(age_kyr, y, t0_grid, tau_rec=TAU_REC, tau_mix=TAU_MIX, min_pre=5, min_post=8): """Matched-filter amplitude at each onset age t0 (kyr). Returns a_hat, formal sigma, n_pre, n_post (samples in first 100 kyr after onset) and the unit response function's self-recovery.""" a = np.full(len(t0_grid), np.nan); se = np.full(len(t0_grid), np.nan) npre = np.zeros(len(t0_grid), int); npost = np.zeros(len(t0_grid), int) for j, t0 in enumerate(t0_grid): lo = np.searchsorted(age_kyr, t0 - POST, "left"); hi = np.searchsorted(age_kyr, t0 + PRE, "right") if hi - lo < min_pre + min_post: continue t = age_kyr[lo:hi]; x = t0 - t # elapsed time since onset (kyr) npre[j] = int(np.sum(x < 0)); npost[j] = int(np.sum((x > 0) & (x <= 100))) if npre[j] < min_pre or npost[j] < min_post: continue g = template(x, tau_rec, tau_mix) X = np.column_stack([np.ones_like(x), x / 100.0, g]) XtX = X.T @ X try: cov = np.linalg.inv(XtX) except np.linalg.LinAlgError: continue beta = cov @ (X.T @ y[lo:hi]) res = y[lo:hi] - X @ beta s2 = res @ res / max(len(x) - 3, 1) a[j] = beta[2]; se[j] = np.sqrt(s2 * cov[2, 2]) return a, se, npre, npost def local_sigma(t0_grid, a, half=2500.0): """Robust (MAD) scatter of the matched-filter amplitude within +-half kyr.""" sig = np.full(len(a), np.nan); med = np.full(len(a), np.nan) ok = np.isfinite(a) tg = t0_grid[ok]; ag = a[ok] for j, t0 in enumerate(t0_grid): lo = np.searchsorted(tg, t0 - half); hi = np.searchsorted(tg, t0 + half) if hi - lo < 50: continue v = ag[lo:hi] m = np.median(v) med[j] = m; sig[j] = 1.4826 * np.median(np.abs(v - m)) return med, sig def background_counts(t0_grid, a, m0=M0, sep_kyr=100.0): """Number of distinct natural negative excursions whose matched-filter amplitude is at least as negative as that of a pulse of M GtC (minima separated by >= sep_kyr).""" rows = [] for M in MASSES: th = float(amp(M, m0)) order = np.argsort(np.where(np.isfinite(a), a, 9.0)) taken = [] for j in order: if not np.isfinite(a[j]) or a[j] > th: break if all(abs(t0_grid[j] - t0_grid[k]) > sep_kyr for k in taken): taken.append(j) ages = sorted(float(t0_grid[j] / 1000.0) for j in taken) rows.append({"pulse_GtC": M, "amplitude_permil": th, "n_natural_excursions_as_large": len(taken), "n_in_0_23Ma": int(sum(x < 23.03 for x in ages)), "n_in_23_66Ma": int(sum(x >= 23.03 for x in ages)), "onset_ages_ma": ";".join(f"{x:.2f}" for x in ages[:60])}) return pd.DataFrame(rows) def run(tau_rec=TAU_REC, tau_mix=TAU_MIX, k_sigma=K_SIGMA, n_inj=4000, seed=SEED, write=True, m0=M0, tag=""): d, n_all = load() age_kyr = d["age_ma"].values * 1000.0 y = d["d13c"].values step = 5.0 t0_grid = np.arange(200.0, 66500.0, step) a, se, npre, npost = scan(age_kyr, y, t0_grid, tau_rec, tau_mix) med, sig = local_sigma(t0_grid, a) z = (a - med) / sig scan_df = pd.DataFrame({"onset_age_ma": t0_grid / 1000.0, "a_hat_permil": a, "a_formal_se": se, "local_median": med, "local_sigma": sig, "z": z, "n_pre": npre, "n_post100": npost}) # ---- natural events: local minima of z below -k --------------------------------- ev = [] zz = np.where(np.isfinite(z), z, 0.0) below = zz < -k_sigma j = 0 while j < len(zz): if below[j]: k = j while k + 1 < len(zz) and (below[k + 1] or (k + 2 < len(zz) and below[k + 2])): k += 1 seg = slice(j, k + 1) jm = j + int(np.argmin(zz[seg])) ev.append({"onset_age_ma": t0_grid[jm] / 1000.0, "a_hat_permil": a[jm], "z": z[jm], "equiv_pulse_GtC": float(mass_from_amp(min(a[jm], -1e-6), m0)), "local_sigma_permil": sig[jm]}) j = k + 1 else: j += 1 # merge events closer than 60 kyr (same excursion picked twice) ev = sorted(ev, key=lambda e: e["onset_age_ma"]) merged = [] for e in ev: if merged and e["onset_age_ma"] - merged[-1]["onset_age_ma"] < 0.06: if e["z"] < merged[-1]["z"]: merged[-1] = e else: merged.append(e) ev_df = pd.DataFrame(merged) # ---- injections ------------------------------------------------------------------- rng = np.random.default_rng(seed) t_inj = rng.uniform(300.0, 66000.0, n_inj) rec = np.zeros((n_inj, len(MASSES)), bool) possible = np.zeros(n_inj, bool) resp_pk = np.full(n_inj, np.nan) for i, ti in enumerate(t_inj): # t0 values tried by the 'observer': within +-15 kyr of the true onset js = np.where(np.abs(t0_grid - ti) <= 15.0)[0] js = js[np.isfinite(a[js]) & np.isfinite(sig[js])] if len(js) == 0: continue possible[i] = True best = np.full(len(MASSES), np.inf) for jj in js: t0 = t0_grid[jj] lo = np.searchsorted(age_kyr, t0 - POST, "left"); hi = np.searchsorted(age_kyr, t0 + PRE, "right") t = age_kyr[lo:hi]; x = t0 - t g = template(x, tau_rec, tau_mix) X = np.column_stack([np.ones_like(x), x / 100.0, g]) sigv = template(ti - t, tau_rec, tau_mix) # injected unit signal at the sample ages beta = np.linalg.lstsq(X, sigv, rcond=None)[0] R = beta[2] # recovered fraction of unit amplitude resp_pk[i] = np.nanmax([resp_pk[i], R]) for m, M in enumerate(MASSES): znew = (a[jj] + amp(M, m0) * R - med[jj]) / sig[jj] best[m] = min(best[m], znew) rec[i] = best < -k_sigma rows = [] for lo, hi, lab in EPOCHS + [(0, 66, "Cenozoic (0.3-66 Ma)"), (0, 23.03, "Neogene+Quaternary")]: sel = (t_inj / 1000.0 >= lo) & (t_inj / 1000.0 < hi) r = {"interval": lab, "age_top_ma": lo, "age_base_ma": hi, "n_injections": int(sel.sum()), "frac_with_adequate_sampling": float(possible[sel].mean()), "median_local_sigma_permil": float(np.nanmedian(sig[(t0_grid / 1000 >= lo) & (t0_grid / 1000 < hi)]))} for m, M in enumerate(MASSES): r[f"P_flag_{M}GtC"] = float(rec[sel, m].mean()) # mass at which half / 90% of injections are flagged (log-interpolated) pm = np.array([rec[sel, m].mean() for m in range(len(MASSES))]) for q, nm in ((0.5, "M50_GtC"), (0.9, "M90_GtC")): if pm.max() < q: r[nm] = np.nan elif pm.min() >= q: r[nm] = MASSES[0] else: r[nm] = float(10 ** np.interp(q, pm, np.log10(MASSES))) rows.append(r) rec_df = pd.DataFrame(rows) # ---- sampling resolution ----------------------------------------------------------- dt = np.diff(age_kyr) mid = 0.5 * (age_kyr[1:] + age_kyr[:-1]) / 1000.0 rr = [] for lo in range(0, 67): s = (mid >= lo) & (mid < lo + 1) if s.sum() == 0: rr.append({"age_ma": lo + 0.5, "n_samples": 0}); continue rr.append({"age_ma": lo + 0.5, "n_samples": int(s.sum()), "median_spacing_kyr": float(np.median(dt[s])), "p90_spacing_kyr": float(np.percentile(dt[s], 90)), "max_gap_kyr": float(dt[s].max())}) res_df = pd.DataFrame(rr) span = age_kyr[-1] - age_kyr[0] summ = { "n_rows_file": int(n_all), "n_d13c": int(len(d)), "age_min_ma": float(d["age_ma"].min()), "age_max_ma": float(d["age_ma"].max()), "n_sites": int(d["site"].nunique()), "median_spacing_kyr": float(np.median(dt)), "mean_spacing_kyr": float(span / len(dt)), "frac_time_in_gaps_gt_10kyr": float(dt[dt > 10].sum() / span), "frac_time_in_gaps_gt_20kyr": float(dt[dt > 20].sum() / span), "frac_time_in_gaps_gt_50kyr": float(dt[dt > 50].sum() / span), "max_gap_kyr": float(dt.max()), "tau_rec_kyr": tau_rec, "tau_mix_kyr": tau_mix, "k_sigma": k_sigma, "M0_GtC": m0, "d13c_source": DSRC, "n_natural_events_flagged": int(len(ev_df)), "amp_permil_by_mass": {str(M): float(amp(M, m0)) for M in MASSES}, "template_peak_fraction": float(template(np.linspace(0, 100, 2001), tau_rec, tau_mix).max()), "median_recovered_fraction_of_injected_amplitude": float(np.nanmedian(resp_pk)), } if write and tag: ev_df.to_csv(os.path.join(DATA, f"cenogrid_events{tag}.csv"), index=False, float_format="%.5g") rec_df.to_csv(os.path.join(DATA, f"cenogrid_recovery{tag}.csv"), index=False, float_format="%.4g") elif write: background_counts(t0_grid, a, m0).to_csv(os.path.join(DATA, "cenogrid_background_counts.csv"), index=False) scan_df.to_csv(os.path.join(DATA, "cenogrid_scan.csv"), index=False, float_format="%.5g") ev_df.to_csv(os.path.join(DATA, "cenogrid_events.csv"), index=False, float_format="%.5g") rec_df.to_csv(os.path.join(DATA, "cenogrid_recovery.csv"), index=False, float_format="%.4g") res_df.to_csv(os.path.join(DATA, "cenogrid_resolution.csv"), index=False, float_format="%.4g") if write and not tag: d[["age_ma", "d13c", "d18o"]].to_csv(os.path.join(DATA, "cenogrid_d13c_series.csv"), index=False, float_format="%.6g") return summ, ev_df, rec_df, res_df STEP_POST = (10.0, 60.0) # kyr after onset (younger side) averaged for the step statistic STEP_PRE = (5.0, 105.0) # kyr before onset (older side) STEP_MASSES = [500, 700, 1000, 1500, 2000, 2500, 3000, 4000, 5000, 10000] def _step_stat(age_kyr, cs, t0, min_post=4, min_pre=6): """Mean d13C in the 50 kyr after an onset at t0 minus the mean in the 100 kyr before it. Ages increase into the past, so 'after' is the younger side.""" a = np.searchsorted(age_kyr, t0 - STEP_POST[1]); b = np.searchsorted(age_kyr, t0 - STEP_POST[0]) c = np.searchsorted(age_kyr, t0 + STEP_PRE[0]); e = np.searchsorted(age_kyr, t0 + STEP_PRE[1]) if b - a < min_post or e - c < min_pre: return np.nan return (cs[b] - cs[a]) / (b - a) - (cs[e] - cs[c]) / (e - c) def run_step(n_inj=1500, seed=SEED + 1, tau_rec=TAU_REC, tau_mix=TAU_MIX, m0=M0): """Injection-recovery with a plain step statistic (no template, no baseline fit). Independent of the matched filter in run(): different statistic, different windows, pulses are added to the data themselves and the statistic recomputed. Used as a robustness check on which natural excursions count as singular and on the detection masses.""" d, _ = load() age_kyr = d["age_ma"].values * 1000.0 y = d["d13c"].values grid = np.arange(300.0, 66500.0, 5.0) cs = np.concatenate([[0.0], np.cumsum(y)]) s0 = np.array([_step_stat(age_kyr, cs, g) for g in grid]) med, sig = local_sigma(grid, s0) z0 = (s0 - med) / sig # response of the statistic to a unit-amplitude template (for 'equivalent pulse' masses) xx = np.arange(-400.0, 400.0, 0.5) # elapsed kyr since onset g_unit = template(xx, tau_rec, tau_mix) resp = g_unit[(xx >= STEP_POST[0]) & (xx < STEP_POST[1])].mean() ev_tabs = {} for k in (5.0, 3.0): idx = np.where(np.isfinite(z0) & (z0 < -k))[0] ev = [] for j in idx: if not ev or grid[j] - grid[ev[-1]] > 150.0: ev.append(j) elif z0[j] < z0[ev[-1]]: ev[-1] = j ev_tabs[k] = pd.DataFrame([{"criterion_sigma": k, "onset_age_ma": grid[j] / 1000.0, "step_permil": s0[j], "z": z0[j], "local_sigma_permil": sig[j], "equiv_pulse_GtC": float(mass_from_amp(min((s0[j] - med[j]) / resp, -1e-6), m0))} for j in ev]) rng = np.random.default_rng(seed) t_inj = rng.uniform(300.0, 66000.0, n_inj) zbest = np.full((n_inj, len(STEP_MASSES)), np.inf) for i, ti in enumerate(t_inj): js = np.where(np.abs(grid - ti) <= 15.0)[0] js = js[np.isfinite(med[js]) & np.isfinite(sig[js])] lo = np.searchsorted(age_kyr, ti - 450.0); hi = np.searchsorted(age_kyr, ti + 150.0) g = template(ti - age_kyr[lo:hi], tau_rec, tau_mix) for m, M in enumerate(STEP_MASSES): y2 = y.copy(); y2[lo:hi] += amp(M, m0) * g c2 = np.concatenate([[0.0], np.cumsum(y2)]) for j in js: st = _step_stat(age_kyr, c2, grid[j]) if np.isfinite(st): zbest[i, m] = min(zbest[i, m], (st - med[j]) / sig[j]) rows = [] for k in (5.0, 3.0): recd = zbest < -k for lo, hi, lab in EPOCHS + [(0, 66, "Cenozoic (0.3-66 Ma)")]: sel = (t_inj / 1000.0 >= lo) & (t_inj / 1000.0 < hi) r = {"criterion_sigma": k, "interval": lab, "n_injections": int(sel.sum())} pm = recd[sel].mean(0) for m, M in enumerate(STEP_MASSES): r[f"P_flag_{M}GtC"] = float(pm[m]) for q, nm in ((0.5, "M50_GtC"), (0.9, "M90_GtC")): r[nm] = float(10 ** np.interp(q, pm, np.log10(STEP_MASSES))) if pm.max() >= q > pm.min() else np.nan rows.append(r) rec = pd.DataFrame(rows) ev = pd.concat([ev_tabs[5.0], ev_tabs[3.0]], ignore_index=True) ev.to_csv(os.path.join(DATA, "cenogrid_step_events.csv"), index=False, float_format="%.5g") rec.to_csv(os.path.join(DATA, "cenogrid_step_recovery.csv"), index=False, float_format="%.4g") return ev, rec, float(np.nanmedian(sig)), float(resp) def spacing_coverage(d): """Share of time lying in gaps between consecutive CENOGRID samples that are no longer than X. This is the sampling density of ONE published splice; it is a floor on what exists.""" age_kyr = d["age_ma"].values * 1000.0 dt = np.diff(age_kyr); mid = 0.5 * (age_kyr[1:] + age_kyr[:-1]) / 1000.0 rows = [] for X in (0.1, 0.3, 1.0, 3.33, 10.0, 50.0): r = {"max_spacing_kyr": X} for lo, hi, lab in ((0.0, 2.58, "0-2.58Ma"), (2.58, 66.0, "2.58-66Ma"), (0.0, 66.0, "0-66Ma")): s = (mid >= lo) & (mid < hi) r[f"covered_myr_{lab}"] = float(dt[s & (dt <= X)].sum() / 1000.0) r[f"frac_time_{lab}"] = r[f"covered_myr_{lab}"] / (hi - lo) r["floor_share_of_2.58-500Ma"] = r["covered_myr_2.58-66Ma"] / (500.0 - 2.58) r["floor_share_of_2.58-2400Ma"] = r["covered_myr_2.58-66Ma"] / (2400.0 - 2.58) rows.append(r) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA, "cenogrid_spacing_coverage.csv"), index=False, float_format="%.4g") return out def main(): summ, ev, rec, res = run() # the looser 'listed among the largest excursions' criterion (3 sigma), full-size run, saved summ3, ev3, rec3, _ = run(k_sigma=3.0, tag="_k3") summ["k3_n_natural_events_flagged"] = int(len(ev3)) summ["k3_events"] = ev3.sort_values("z").to_dict("records") # robustness: other recovery times, mixing times, thresholds, reservoir size rob = [] for lab, kw in (("baseline", {}), ("tau_rec=50kyr", {"tau_rec": 50.0}), ("tau_rec=200kyr", {"tau_rec": 200.0}), ("tau_mix=3kyr", {"tau_mix": 3.0}), ("tau_mix=15kyr", {"tau_mix": 15.0}), ("k=3", {"k_sigma": 3.0}), ("k=4", {"k_sigma": 4.0}), ("k=6", {"k_sigma": 6.0}), ("M0=60000GtC", {"m0": 60000.0})): if lab == "baseline": s2, e2, r2 = summ, ev, rec elif lab == "k=3": s2, e2, r2 = summ3, ev3, rec3 else: s2, e2, r2, _ = run(write=False, n_inj=1500, **kw) row = {"variant": lab, "n_natural_flagged": int(len(e2))} for iv in ("Cenozoic (0.3-66 Ma)", "Neogene+Quaternary", "Eocene", "Paleocene"): rr = r2[r2["interval"] == iv].iloc[0] row[f"M50_{iv}"] = rr["M50_GtC"]; row[f"M90_{iv}"] = rr["M90_GtC"] row[f"P700_{iv}"] = rr["P_flag_700GtC"]; row[f"P5000_{iv}"] = rr["P_flag_5000GtC"] rob.append(row) rob = pd.DataFrame(rob) rob.to_csv(os.path.join(DATA, "cenogrid_robustness.csv"), index=False, float_format="%.4g") summ["events_top"] = ev.sort_values("z").head(15).to_dict("records") # ---- alternative detector and the record's own sampling density ----------------------- sev, srec, ssig, sresp = run_step() summ["step_detector"] = { "windows_kyr": {"post": STEP_POST, "pre": STEP_PRE}, "n_injections": 1500, "median_local_sigma_permil": ssig, "response_to_unit_template": sresp, "events_5sigma": sev[sev["criterion_sigma"] == 5.0].to_dict("records"), "events_3sigma": sev[sev["criterion_sigma"] == 3.0].to_dict("records"), "cenozoic": {f"{int(k)}sigma": srec[(srec["criterion_sigma"] == k) & (srec["interval"] == "Cenozoic (0.3-66 Ma)")].iloc[0].to_dict() for k in (5.0, 3.0)}} # largest natural excursions other than the PETM, by equivalent pulse mass, matched filter (3-sigma list) others = ev3[(ev3["onset_age_ma"] < 55.5) | (ev3["onset_age_ma"] > 56.5)].sort_values("equiv_pulse_GtC", ascending=False) summ["largest_non_PETM_equiv_pulse_GtC_matched_filter"] = others[["onset_age_ma", "equiv_pulse_GtC", "z"]].head(4).to_dict("records") d_all, _ = load() sc = spacing_coverage(d_all) summ["spacing_coverage"] = sc.to_dict("records") summ["samples_per_myr_0_66Ma"] = float(((d_all["age_ma"] < 66)).sum() / 66.0) r5 = rec[rec["interval"] == "Cenozoic (0.3-66 Ma)"].iloc[0] summ["recorded_peak_permil_by_mass"] = {str(M): float(amp(M)) * summ["template_peak_fraction"] for M in MASSES} summ["P_flag_3000GtC_by_epoch_5sigma"] = {r["interval"]: float(r["P_flag_3000GtC"]) for _, r in rec.iterrows()} summ["P_flag_500GtC_cenozoic"] = {"5sigma": float(r5["P_flag_500GtC"]), "3sigma": float(rec3[rec3["interval"] == "Cenozoic (0.3-66 Ma)"].iloc[0]["P_flag_500GtC"])} save_json(summ, "cenogrid_summary.json") print(sev.to_string()); print(srec.to_string()); print(sc.to_string()) pd.set_option("display.width", 250) print(json.dumps({k: v for k, v in summ.items() if k not in ("events_top", "k3_events")}, indent=1)) print(pd.read_csv(os.path.join(DATA, "cenogrid_background_counts.csv")).iloc[:, :5].to_string()) print(ev3.sort_values("z").to_string()); print(rec3.to_string()) print(ev.sort_values("z").head(40).to_string()) print(rec.to_string()) print(rob.to_string()) if __name__ == "__main__": main()