"""The controls and the floor, run before the real statistic is computed. Nothing here prints the syzygy excess of the real catalogue. It prints the time measure, the two null controls, the declustered size, and the floor -- the smallest injected excess this instrument can recover. If the floor is above the +5 % bar, the SURVIVED branch of the sealed kill rule cannot fire on any data, and that is written down here rather than discovered afterwards. python controls.py > controls.txt The phase grid is cached in ~/data/quake-moon/ because pyephem is one call per instant and the fine grid is four million of them. """ import csv import pathlib import numpy as np import rule HERE = pathlib.Path(__file__).parent DATA = pathlib.Path.home() / "data" / "quake-moon" CAT = DATA / "comcat_m60_1973_2025.csv" rng = np.random.default_rng(rule.SEED) T0 = np.datetime64(rule.START[:-1], "s") T1 = np.datetime64(rule.END[:-1], "s") SPAN = int((T1 - T0) / np.timedelta64(1, "s")) def load(): times, lat, lon, dep, mag = [], [], [], [], [] with CAT.open() as f: for row in csv.DictReader(f): if not row["time"] or not row["mag"]: continue times.append(np.datetime64(row["time"][:19], "s")) lat.append(float(row["latitude"])) lon.append(float(row["longitude"])) dep.append(float(row["depth"]) if row["depth"] else np.nan) mag.append(float(row["mag"])) return (np.array(times), np.array(lat), np.array(lon), np.array(dep), np.array(mag)) def grid_phases(step_seconds, tag): cache = DATA / f"grid_{tag}.npy" if cache.exists(): return np.load(cache) n = SPAN // step_seconds + 1 g = T0 + (np.arange(n) * step_seconds).astype("timedelta64[s]") th = rule.phases(g) np.save(cache, th) return th def ci_normal(inside, f_exp): """The normal approximation to the row bootstrap, used where 10,000 percentile draws x 200 simulations would not finish. Checked against the real percentile bootstrap in analyse.py.""" n = inside.size f = inside.mean() se = np.sqrt(f * (1 - f) / n) / f_exp e = f / f_exp - 1.0 return e, e - 1.96 * se, e + 1.96 * se def main(): times, lat, lon, dep, mag = load() print(f"catalogue: {len(times)} events, M>=6.0, " f"{str(times.min())[:10]} .. {str(times.max())[:10]}") # ---- time-measure control th_hour = grid_phases(3600, "hour") th_fine = grid_phases(420, "7min") f_hour = float((rule.syzygy_distance(th_hour) <= rule.HALF).mean()) f_fine = float((rule.syzygy_distance(th_fine) <= rule.HALF).mean()) print(f"\ntime measure hourly grid f_exp = {f_hour:.6f} " f"({th_hour.size} points)") print(f"time measure 7-min grid f_exp = {f_fine:.6f} " f"({th_fine.size} points)") print(f" agreement |hour - fine| = {abs(f_hour - f_fine):.6f} " f"(rule: < 0.001) -> {'PASS' if abs(f_hour-f_fine) < 0.001 else 'FAIL'}") print(f" departure from 0.5 = {abs(f_hour - 0.5):.6f} " f"(rule: < 0.02) -> {'PASS' if abs(f_hour-0.5) < 0.02 else 'FAIL'}") f_exp = f_hour f_full = float((rule.full_distance(th_hour) <= rule.HALF).mean()) print(f"time measure full moon alone (S3): f_exp = {f_full:.6f}") # ---- declustering mains, cluster = rule.decluster(times, lat, lon, mag) print(f"\ndeclustering (Gardner-Knopoff, forward windows)") print(f" events {len(times)}") print(f" mainshocks {len(mains)} ({100*len(mains)/len(times):.1f} %)") print(f" clusters {len(np.unique(cluster))}") big = np.bincount(cluster).max() print(f" largest cluster {big} events") print(f" rule: fewer than 1000 mainshocks voids the study -> " f"{'PASS' if len(mains) >= 1000 else 'FAIL'}") n_main = len(mains) mains_s, _ = rule.decluster(times, lat, lon, mag, symmetric=True) print(f" symmetric variant (S4): {len(mains_s)} mainshocks") # ---- null control N1: uniform times through the real phase code print(f"\nN1 times drawn uniformly from the window, N = {n_main}") draws = 200 es, excl = [], 0 for _ in range(draws): t = T0 + rng.integers(0, SPAN, size=n_main).astype("timedelta64[s]") th = rule.phases(t) inside = (rule.syzygy_distance(th) <= rule.HALF).astype(float) e, lo, hi = ci_normal(inside, f_exp) es.append(e) excl += (lo > 0) or (hi < 0) es = np.array(es) print(f" mean E over {draws} draws = {100*es.mean():+.3f} % " f"(sd {100*es.std():.3f})") print(f" 'interval excludes 0' fired {excl}/{draws} = " f"{100*excl/draws:.1f} % (rule: near 5 %)") print(f" -> {'PASS' if abs(es.mean()) < 0.01 and excl/draws < 0.12 else 'FAIL'}") # ---- null control N2: jitter the real times by +/- 365 days print(f"\nN2 every real mainshock time jittered by U(-365, +365) days") tm = times[mains] es2 = [] for _ in range(50): j = rng.integers(-365 * 86400, 365 * 86400, size=n_main) t = tm + j.astype("timedelta64[s]") th = rule.phases(t) inside = (rule.syzygy_distance(th) <= rule.HALF).astype(float) es2.append(ci_normal(inside, f_exp)[0]) es2 = np.array(es2) print(f" mean E over 50 draws = {100*es2.mean():+.3f} % " f"(sd {100*es2.std():.3f})") print(f" -> {'PASS' if abs(es2.mean()) < 0.01 else 'FAIL'}") # ---- the floor print(f"\nFLOOR injected excess recovered, {n_main} mainshocks, " f"200 simulations each") floor = None for delta in (0.02, 0.05, 0.10, 0.20): p = f_exp * (1 + delta) hits = 0 for _ in range(200): inside = (rng.random(n_main) < p).astype(float) e, lo, hi = ci_normal(inside, f_exp) hits += lo > 0 rate = hits / 200 mark = "" if floor is None and rate >= 0.90: floor = delta mark = " <- floor" print(f" +{100*delta:>5.1f} % recovered (interval excludes 0) " f"{100*rate:5.1f} % of the time{mark}") print(f" floor = {'none of those' if floor is None else f'+{100*floor:.0f} %'}" f" bar = +{100*rule.BAR:.0f} %") if floor is None or floor > rule.BAR: print(" -> the SURVIVED branch CANNOT fire: void in that direction") else: print(" -> the floor is at or below the bar: both branches live") # ---- the same, for secondary S2 (M >= 7.0), before its data is read n7 = int((mag[mains] >= 7.0).sum()) print(f"\nFLOOR for S2 (M>=7.0): {n7} mainshocks") f7 = None for delta in (0.05, 0.10, 0.20, 0.40): p = f_exp * (1 + delta) hits = sum(ci_normal((rng.random(n7) < p).astype(float), f_exp)[1] > 0 for _ in range(200)) rate = hits / 200 if f7 is None and rate >= 0.90: f7 = delta print(f" +{100*delta:>5.1f} % recovered {100*rate:5.1f} %") print(f" floor = {'none of those' if f7 is None else f'+{100*f7:.0f} %'}" f" -> S2 is {'VOID: unfailable' if f7 is None or f7 > rule.BAR else 'live'}") np.save(DATA / "f_exp.npy", np.array([f_exp, f_full])) if __name__ == "__main__": main()