"""Run the sealed rule on the real series.

    python analyse.py > results.txt

Needs ~/data/births/ from fetch.py. The SSA file 2000-2014 carries the
verdict; the NCHS file, restricted to 1994-1999, is a second witness and is
printed as a description only. Writes curve.csv (mean residual by day of the
lunar month, both series) beside this script.
"""

import datetime as dt
import json
import pathlib

import pandas as pd

import rule

DATA = pathlib.Path.home() / "data" / "births"
HERE = pathlib.Path(__file__).resolve().parent


def _plain(v):
    if isinstance(v, (tuple, list)):
        return [_plain(x) for x in v]
    if isinstance(v, (bool, str, int)):
        return v
    return round(float(v), 4)


def load(name, first=None, last=None):
    df = pd.read_csv(DATA / name)
    df["date"] = [dt.date(int(y), int(m), int(d))
                  for y, m, d in zip(df["year"], df["month"], df["date_of_month"])]
    if first:
        df = df[df["year"] >= first]
    if last:
        df = df[df["year"] <= last]
    return df[["date", "births"]]


def report(label, out):
    print(f"== {label}")
    print(f"days {out['days']}, kept in the fit {out['kept']}, "
          f"window days kept {out['window_days']}, residual SD {100*out['residual_sd']:.2f} %")
    print(f"full-moon date + day after: {out['est_pct']:+.3f} % "
          f"(95 % interval {out['lo_pct']:+.3f} % to {out['hi_pct']:+.3f} %), "
          f"SE {100*out['se']:.3f} % (log scale)")
    print(f"floor (smallest excess called survived 4 times in 5): {out['floor_pct']:.3f} %")
    v, h = out["valentine"], out["halloween"]
    print(f"control, 14 Feb vs a week either side: {v[0]:+.2f} % ({v[1]:+.2f} to {v[2]:+.2f}), {v[3]} years")
    print(f"control, 31 Oct vs a week either side: {h[0]:+.2f} % ({h[1]:+.2f} to {h[2]:+.2f}), {h[3]} years")
    print(f"controls pass: {out['controls_pass']}")
    print(f"verdict under the sealed rule: {out['verdict']}")
    print()


def main():
    ssa = rule.analyse(load("US_births_2000-2014_SSA.csv"))
    report("SSA 2000-2014 (sealed)", ssa)
    nchs = rule.analyse(load("US_births_1994-2003_CDC_NCHS.csv", 1994, 1999))
    report("CDC/NCHS 1994-1999 (second witness, described)", nchs)
    curve = pd.DataFrame({"day_of_lunar_month": sorted(ssa["curve_pct"]),
                          "ssa_2000_2014_pct": [ssa["curve_pct"][k] for k in sorted(ssa["curve_pct"])]})
    curve["nchs_1994_1999_pct"] = curve["day_of_lunar_month"].map(nchs["curve_pct"])
    curve.round(4).to_csv(HERE / "curve.csv", index=False)
    keep = ("days", "kept", "window_days", "residual_sd", "est_pct", "lo_pct", "hi_pct",
            "floor_pct", "verdict", "valentine", "halloween", "controls_pass")
    json.dump({"ssa_2000_2014": {k: _plain(ssa[k]) for k in keep},
               "nchs_1994_1999": {k: _plain(nchs[k]) for k in keep}},
              (HERE / "results.json").open("w"), indent=1)
    print("mean residual by day of the lunar month (0 = full-moon date), SSA:")
    print(" ".join(f"{k}:{v:+.2f}" for k, v in sorted(ssa["curve_pct"].items())))


if __name__ == "__main__":
    main()