"""The 2000-2009 and 2013-2025 blocks taken together (exploratory; the pre-registered test scores each on its own).

Each block's exact null distribution comes from overlap_test.py's circular-shift test (every shift of the block's dates
once); the blocks are shifted independently, so the null distribution of their summed overlaps is the convolution
of the two. Also gives the ratio of observed to expected overlaps with an approximate 95% interval (exact Poisson,
from intervals.py). Standard library only.

Usage: python3 -I combine_blocks.py <dates json A> <usgs csv A> <dates json B> <usgs csv B> [--out file.json]
"""
import datetime, json, pathlib, sys
from collections import Counter
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import overlap_test as ot
import intervals as iv


def block(dates_path, csv_path):
    dates = json.load(open(dates_path, encoding="utf-8"))
    years = sorted(int(y) for y in dates)
    quakes = ot.quake_days([csv_path], set(years))
    start = datetime.date(years[0], 1, 1); L = (datetime.date(years[-1] + 1, 1, 1) - start).days
    A = {ot.index(d, start) for y in years for d in dates[str(y)]}
    B = {ot.index(d, start) for d in quakes}
    dist = [sum(((a + s - 1) % L) + 1 in B for a in A) for s in range(L)]
    return {"years": f"{years[0]}-{years[-1]}", "L": L, "dates": len(A), "quake_days": len(B),
            "observed": dist[0], "expected": sum(dist) / L}, dist


def main():
    args = sys.argv[1:]; out = None
    if "--out" in args:
        i = args.index("--out"); out = args[i + 1]; del args[i:i + 2]
    (a, da), (b, db) = block(args[0], args[1]), block(args[2], args[3])
    ca, cb = Counter(da), Counter(db)
    k = a["observed"] + b["observed"]; e = a["expected"] + b["expected"]
    p = sum(na * nb for ka, na in ca.items() for kb, nb in cb.items() if ka + kb >= k) / (a["L"] * b["L"])
    lo, hi = iv.poisson_ci(k)
    res = {"blocks": [a, b], "observed": k, "expected": round(e, 3), "ratio": round(k / e, 3),
           "ratio_95_interval_approx": [round(lo / e, 2), round(hi / e, 2)], "p_exact_convolved": round(p, 5)}
    for x in res["blocks"]:
        x["expected"] = round(x["expected"], 3)
    print(json.dumps(res, indent=1))
    if out:
        with open(out, "w", encoding="utf-8", newline="\n") as f:
            json.dump(res, f, indent=1); f.write("\n")


if __name__ == "__main__":
    main()