"""Collect the numbers on /planet-alignments-earthquakes/ into one small JSON file (standard library only). Reads the outputs of select_dates.py, overlap_test.py, outsample_test.py, diagnose_dates.py, explore.py, why_2011.py, intervals.py and tides.py from , and the Horizons tables for the SHA-256 of their ephemeris lines (the text between $$SOE and $$EOE, which is the same on every fetch; the header carries the time of the request, so whole-file sums change). Usage: python3 -I summarize.py """ import hashlib, json, pathlib, sys FROZEN = ["fetch_horizons.py", "fetch_usgs.py", "config_dates.py", "select_dates.py", "overlap_test.py", "outsample_test.py", "diagnose_dates.py"] VARIANTS = [("", "The paper's choices: the Sun's row first, the SVD, alpha = 0.35"), ("-a0.35-svd-sunlast", "The Sun's row listed last"), ("-a0.35-svd-reversed", "The bodies listed in reverse order after the Sun"), ("-a0.35-eigh-sunfirst", "NumPy's symmetric eigensolver (eigh) in place of the SVD"), ("-a1.0-svd-sunfirst", "The threshold eps alone (alpha = 1)")] def load(p): return json.load(open(p, encoding="utf-8")) def main(): rdir, hdir, cdir, out = map(pathlib.Path, sys.argv[1:5]) res = {"about": "Numbers behind https://factchecked.science/planet-alignments-earthquakes/, from our rebuild of " "Nikitopoulos (2026) on JPL Horizons and USGS ComCat data fetched on 2026-10-28. Dates are TDB " "calendar dates; earthquake days are UTC days with a USGS event of magnitude above 6.5."} sel = {} for tag in ("2000-2009", "2010-2012", "2013-2025"): for yr in load(rdir / f"selection-{tag}.json"): sel[yr["year"]] = yr ins = load(rdir / "overlap-2010-2012.json") blocks = {"2010-2012": ins["2010-2012"], "2013-2025": load(rdir / "outsample-2013-2025.json"), "2000-2009": load(rdir / "outsample-2000-2009.json")} rows = [] for y in range(2000, 2026): if 2010 <= y <= 2012: r, block = ins[str(y)], "2010-2012 (the paper's years)" else: tag = "2013-2025" if y >= 2013 else "2000-2009" r, block = blocks[tag]["per_year"][str(y)], tag + (" (primary test)" if y >= 2013 else " (second block)") s = sel[y] rows.append({"year": y, "block": block, "bodies": ["Sun"] + s["selected"], "candidates_by_step": [st["candidates"] for st in s["steps"]], "dates": r["m_dates"], "quake_days": r["n_quake_days"], "overlap": r["observed"], "expected": r["expected_exact"], "p_exact": r["p_exact"]}) res["per_year"] = rows keep = ("L", "m_dates", "n_quake_days", "observed", "expected_exact", "p_exact", "shifts_at_or_above_observed", "q95_exact", "k_crit_p_below_0.05", "power_at_in_sample_rate", "power_at_1.5x_chance", "years_p_below_0.05", "ratio_obs_exp") res["blocks"] = {k: {kk: v[kk] for kk in keep if kk in v} for k, v in blocks.items()} res["blocks"]["2010-2012"]["ratio_obs_exp"] = round(ins["2010-2012"]["observed"] / ins["2010-2012"]["expected_exact"], 3) res["blocks"]["2010-2012"]["p_monte_carlo_200000"] = ins["2010-2012"]["p_mc"] for y in ("2010", "2011", "2012"): res["blocks"].setdefault("2010-2012_by_year", {})[y] = {"p_monte_carlo_200000": ins[y]["p_mc"], "q95_exact": ins[y]["q95_exact"], "overlap_days": ins[y]["overlap_days"]} var = [] paper_dates = {d for ds in load(rdir / "dates-2010-2012.json").values() for d in ds} for suffix, label in VARIANTS: o = load(rdir / f"overlap-2010-2012{suffix}.json") vd = {d for ds in load(rdir / f"dates-2010-2012{suffix}.json").values() for d in ds} var.append({"variant": label, "dates_by_year": [o[y]["m_dates"] for y in ("2010", "2011", "2012")], "dates_shared_with_paper": len(vd & paper_dates), "dates": o["2010-2012"]["m_dates"], "overlap": o["2010-2012"]["observed"], "expected": o["2010-2012"]["expected_exact"], "p_exact": o["2010-2012"]["p_exact"], "p_by_year": [o[y]["p_exact"] for y in ("2010", "2011", "2012")]}) res["variants_2010_2012"] = var flips, align = {}, {} for tag in ("2010-2012", "2013-2025", "2000-2009"): dg = load(rdir / f"diagnose-{tag}.json") resid = [o["after_undoing"]["residual_AU"] for o in dg]; jumps = [o["jump_AU"] for o in dg] flips[tag] = {"candidate_minutes": len(dg), "pure_sign_flips": sum(x < 3.5e-5 for x in resid), "axes_flipped_1": sum(sum(s < 0 for s in o["after_undoing"]["signs"]) == 1 for o in dg), "axes_flipped_2": sum(sum(s < 0 for s in o["after_undoing"]["signs"]) == 2 for o in dg), "reordered": sum(list(o["after_undoing"]["perm"]) != [0, 1, 2] for o in dg), "residual_AU_range": [float(f"{min(resid):.2g}"), float(f"{max(resid):.2g}")], "jump_AU_range": [float(f"{min(jumps):.2g}"), float(f"{max(jumps):.2g}")]} by_date = {} for o in dg: by_date.setdefault(o["stamp"][:10], o["R2_pct_of_year"]) v = list(by_date.values()) align[tag] = {"dates": len(v), "mean_rank_R2_on_a_line": round(sum(v) / len(v), 3), "se_if_unrelated": round((1 / 12 / len(v)) ** 0.5, 3)} for b in load(rdir / "explore.json")["blocks"]: align[b["block"]]["mean_rank_R1_same_side"] = b["mean_rank_R1"] res["sign_flips"] = flips res["alignment_rank"] = align res["lunar_phase"] = load(rdir / "explore.json") res["why_2011"] = load(rdir / "why-2011.json") res["ratio_intervals_exploratory"] = load(rdir / "intervals.json") res["tides"] = load(rdir / "tides.json") res["interpolation_check"] = "Cubic Hermite positions from hourly states agree with Horizons' own 1-minute " \ "output to within 2.8e-13 AU on 2010-01-06, 2010-04-30, 2011-09-16 and 2012-10-28, " \ "with the same candidate minutes (check_interpolation.py)." res["dates"] = {str(y): sel[y]["dates"] for y in range(2000, 2026)} res["combined_23_years_exploratory"] = {"preprint_text_rule": load(rdir / "combined-preregistered-rule.json")} ar = rdir / "author-rule" if ar.exists(): # the further, exploratory test with the author's code's rule asel = {} for tag in ("2000-2009", "2010-2012", "2013-2025"): for yr in load(ar / f"selection-{tag}.json"): asel[yr["year"]] = yr ains = load(ar / "overlap-2010-2012.json") ablk = {"2013-2025": load(ar / "outsample-2013-2025.json"), "2000-2009": load(ar / "outsample-2000-2009.json")} arows = [] for y in range(2000, 2026): r = ains[str(y)] if 2010 <= y <= 2012 else ablk["2013-2025" if y >= 2013 else "2000-2009"]["per_year"][str(y)] arows.append({"year": y, "bodies": ["Sun"] + asel[y]["selected"], "candidates_by_step": [st["candidates"] for st in asel[y]["steps"]], "dates": r["m_dates"], "quake_days": r["n_quake_days"], "overlap": r["observed"], "expected": r["expected_exact"], "p_exact": r["p_exact"], "body_set_differs_from_text_rule": asel[y]["selected"] != sel[y]["selected"]}) aflips, aalign = {}, {} for tag in ("2010-2012", "2013-2025", "2000-2009"): dg = load(ar / f"diagnose-{tag}.json") resid = [o["after_undoing"]["residual_AU"] for o in dg] aflips[tag] = {"candidate_minutes": len(dg), "pure_sign_flips": sum(x < 3.5e-5 for x in resid), "with_a_sign_reversed": sum(any(s < 0 for s in o["after_undoing"]["signs"]) for o in dg), "not_pure": [{"stamp": o["stamp"], "residual_AU": float(f'{o["after_undoing"]["residual_AU"]:.2g}'), "lambda2_over_lambda1": round(o["eigenvalues_AU2"][1] / o["eigenvalues_AU2"][0], 2)} for o in dg if o["after_undoing"]["residual_AU"] >= 3.5e-5], "reordered": sum(list(o["after_undoing"]["perm"]) != [0, 1, 2] for o in dg), "residual_AU_range": [float(f"{min(resid):.2g}"), float(f"{max(resid):.2g}")]} by_date = {} for o in dg: by_date.setdefault(o["stamp"][:10], o["R2_pct_of_year"]) v = list(by_date.values()) aalign[tag] = {"dates": len(v), "mean_rank_R2_on_a_line": round(sum(v) / len(v), 3), "se_if_unrelated": round((1 / 12 / len(v)) ** 0.5, 3)} avar = [] for method, order, label in (("svd", "sunlast", "The Sun's row listed last"), ("svd", "reversed", "The bodies listed in reverse order after the Sun"), ("eigh", "sunfirst", "NumPy's symmetric eigensolver (eigh) in place of the SVD")): vd = ar / f"{method}-{order}" o = load(vd / "overlap-2010-2012.json") entry = {"variant": label, "2010-2012": {k: o["2010-2012"][k] for k in ("m_dates", "observed", "expected_exact", "p_exact")}} for tag in ("2013-2025", "2000-2009"): b = load(vd / f"outsample-{tag}.json") entry[tag] = {k: b[k] for k in ("m_dates", "observed", "expected_exact", "p_exact", "ratio_obs_exp")} avar.append(entry) res["author_rule_exploratory"] = { "rule": "The author's planetary_dates.m (Zenodo 10.5281/zenodo.21737429, v1.1) tests the threshold on Earth's " "row of the embedding only (line 222) and never tests a year's first minute; the preprint's " "section 2.2, and the pre-registered test, compare every coordinate. Run after the pre-registered test.", "per_year": arows, "blocks": {k: {kk: v[kk] for kk in keep if kk in v} for k, v in ablk.items()}, "combined_23_years": load(ar / "combined.json"), "ratio_intervals": load(ar / "intervals.json"), "sign_flips": aflips, "alignment_rank": aalign, "variants": avar, "dates": {str(y): asel[y]["dates"] for y in range(2000, 2026)}} res["frozen_code_sha256"] = {f: hashlib.sha256((cdir / f).read_bytes()).hexdigest() for f in FROZEN} hz = {} for p in sorted(hdir.glob("*_1h.txt")): t = p.read_text(encoding="utf-8") hz[p.name] = hashlib.sha256(t.split("$$SOE", 1)[1].split("$$EOE", 1)[0].encode("utf-8")).hexdigest() res["horizons_ephemeris_lines_sha256"] = hz with open(out, "w", encoding="utf-8", newline="\n") as f: json.dump(res, f, indent=1, ensure_ascii=False); f.write("\n") print("wrote", out, out.stat().st_size, "bytes;", len(hz), "Horizons tables") if __name__ == "__main__": main()