"""Why the selection rule added the Moon and Jupiter in 2011 (exploratory; not part of the frozen test).

With the Sun and the four inner planets, 2011 gives 3,354 candidate minutes. Most of them form one unbroken run,
while the ratio of the Gram matrix's second to first eigenvalue (how evenly the five bodies are spread in the plane:
near 0 for bodies strung along a line, near 1 for an even spread) is near its highest of the year. The script prints
the runs, the ratio on the day the long run starts against its distribution over the year (one value a day, at
00:00 TDB), and the count of candidates with the threshold eps alone (alpha = 1).

Needs NumPy. Usage: python3 -I why_2011.py <horizons dir> [out.json]
"""
import datetime, json, pathlib, sys
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
import numpy as np
import config_dates as cd

INNER = ["199", "299", "399", "499"]


def main():
    hdir = sys.argv[1]; year = 2011
    stamps, jumps = cd.candidates(hdir, year, INNER)
    k = np.nonzero(jumps >= 0.35 * cd.EPS)[0]
    runs = np.split(k, np.nonzero(np.diff(k) > 1)[0] + 1)
    longest = max(runs, key=len)
    pos = [cd.minute_positions(hdir, b, year) for b in INNER]
    X = np.stack([np.zeros_like(pos[0])] + pos, axis=1)
    _, S = cd.embed(X[1::1440])                     # 00:00 TDB of every day
    ratio = S[:, 1] / S[:, 0]
    day = int(longest[0]) // 1440
    start = datetime.datetime(year, 1, 1) + datetime.timedelta(minutes=int(longest[0]))
    alpha1 = int((jumps >= cd.EPS).sum())
    out = {"year": year, "bodies": "Sun, Mercury, Venus, Earth, Mars", "candidates": int(len(k)),
           "runs": len(runs), "longest_run_minutes": int(len(longest)),
           "longest_run_start_TDB": start.strftime("%Y-%m-%d %H:%M"),
           "candidates_outside_longest_run": int(len(k) - len(longest)),
           "lambda2_over_lambda1_on_that_day": round(float(ratio[day]), 3),
           "lambda2_over_lambda1_median_of_year": round(float(np.median(ratio)), 3),
           "lambda2_over_lambda1_99th_percentile": round(float(np.percentile(ratio, 99)), 3),
           "rank_of_that_day_in_year": round(float((ratio < ratio[day]).mean()), 3),
           "candidates_with_alpha_1": alpha1}
    print(json.dumps(out, indent=1))
    if len(sys.argv) > 2:
        with open(sys.argv[2], "w", encoding="utf-8", newline="\n") as f:
            json.dump(out, f, indent=1); f.write("\n")


if __name__ == "__main__":
    main()