"""Rebuild Nikitopoulos (2026) "planetary configuration dates" from JPL Horizons vectors.

Method, from the preprint's section 2.2 (pp. 4-5):
  * heliocentric ecliptic (ICRF) positions in AU at 1-minute steps; the Sun at (0,0,0)
    and the selected bodies form X(t) (N x 3);
  * D_ij = |x_i - x_j|^2; row means Dbar_i; Dhat = (1/N^2) sum_{j>k} D_jk;
    Do_i = Dbar_i - Dhat; G_ij = (Do_i + Do_j - D_ij)/2 (classical MDS Gram matrix);
  * SVD G = U S U^T; X* = U_3 S_3^(1/2);
  * a minute t is a candidate if max_ij |X*_ij(t) - X*_ij(t-1)| >= alpha*eps,
    eps = 1e-4 AU, alpha = 0.35.
Positions at each minute come from cubic Hermite interpolation of Horizons' hourly
positions and velocities (checked against Horizons' own 1-minute output: see
check_interpolation.py). Needs NumPy (the SVD is LAPACK's, as in MATLAB).

Usage: python -I config_dates.py <horizons dir> <year> <body,body,...> [alpha] [svd|eigh] [order]
  bodies in matrix order after the Sun, e.g. 199,299,399,499
  order: 'sunfirst' (default), 'sunlast' (the Sun's row last) or 'reversed' (Sun first,
  then the bodies in reverse order). Only the order of the rows changes, not the geometry.
"""
import calendar, datetime, json, pathlib, sys
import numpy as np

EPS = 1e-4
NAMES = {"199": "Mercury", "299": "Venus", "399": "Earth", "301": "Moon", "499": "Mars",
         "599": "Jupiter", "699": "Saturn", "799": "Uranus", "899": "Neptune", "3": "EMB"}


def load_states(path):
    """Return (jd, states[n,6]) from a Horizons CSV vector table (VEC_TABLE=2)."""
    text = pathlib.Path(path).read_text(encoding="utf-8")
    body = text.split("$$SOE", 1)[1].split("$$EOE", 1)[0]
    jd, st = [], []
    for line in body.strip().splitlines():
        f = [x.strip() for x in line.split(",")]
        jd.append(float(f[0])); st.append([float(v) for v in f[2:8]])
    return np.array(jd), np.array(st)


def minute_positions(hdir, body, year):
    """Positions at every minute of `year` (TDB), with one extra minute before Jan 1 00:00.
    Returns array [nmin+1, 3]; row 0 is Dec 31 23:59 of the previous year."""
    path = pathlib.Path(hdir) / f"{body}_{year - 1}-12-31_{year + 1}-01-01_1h.txt"
    jd, st = load_states(path)
    jd0 = 2451544.5 + (datetime.date(year, 1, 1) - datetime.date(2000, 1, 1)).days  # Jan 1 00:00
    assert abs(jd[0] - (jd0 - 1 / 24)) < 1e-6 and np.allclose(np.diff(jd), 1 / 24, atol=1e-7)
    nmin = (366 if calendar.isleap(year) else 365) * 1440
    k = np.arange(-1, nmin)              # minute index from Jan 1 00:00
    j = (k + 60) // 60                   # hourly knot at or before the minute
    s = ((k + 60) % 60) / 60.0
    h = 1 / 24
    p0, v0, p1, v1 = st[j, :3], st[j, 3:], st[j + 1, :3], st[j + 1, 3:]
    s = s[:, None]
    h00 = 2 * s**3 - 3 * s**2 + 1; h10 = s**3 - 2 * s**2 + s
    h01 = -2 * s**3 + 3 * s**2;    h11 = s**3 - s**2
    return h00 * p0 + h10 * h * v0 + h01 * p1 + h11 * h * v1


def embed(X, method="svd"):
    """X: [T, N, 3] -> X*: [T, N, 3] via the paper's Gram matrix and SVD."""
    diff = X[:, :, None, :] - X[:, None, :, :]
    D = np.einsum("tijk,tijk->tij", diff, diff)
    N = X.shape[1]
    Dbar = D.mean(axis=2)
    il = np.tril_indices(N, -1)
    Dhat = D[:, il[0], il[1]].sum(axis=1) / N**2
    Do = Dbar - Dhat[:, None]
    G = 0.5 * (Do[:, :, None] + Do[:, None, :] - D)
    if method == "svd":
        U, S, _ = np.linalg.svd(G)
        return U[:, :, :3] * np.sqrt(S[:, None, :3]), S
    w, V = np.linalg.eigh(G)             # ascending eigenvalues
    V = V[:, :, ::-1]; w = w[:, ::-1]
    return V[:, :, :3] * np.sqrt(np.clip(w[:, None, :3], 0, None)), w


def candidates(hdir, year, bodies, alpha=0.35, method="svd", order="sunfirst", chunk=44640):
    pos = [minute_positions(hdir, b, year) for b in bodies]
    T = pos[0].shape[0]
    sun = np.zeros_like(pos[0])
    cols = {"sunfirst": [sun] + pos, "sunlast": pos + [sun], "reversed": [sun] + pos[::-1]}[order]
    X = np.stack(cols, axis=1)           # [T, N, 3]
    thr = alpha * EPS
    flagged, jumps = [], np.empty(T - 1)
    prev = None
    for a in range(0, T, chunk):
        b = min(T, a + chunk + 1)
        Xs, _ = embed(X[a:b], method)
        d = np.abs(np.diff(Xs, axis=0)).max(axis=(1, 2))   # jump from minute i-1 to i
        jumps[a:b - 1] = d
    k = np.nonzero(jumps >= thr)[0]       # jumps[k] is between grid rows k and k+1
    jan1 = datetime.datetime(year, 1, 1)
    stamps = [jan1 + datetime.timedelta(minutes=int(i)) for i in k]   # row k+1 = minute k
    return stamps, jumps


def main():
    hdir, year, bodies = sys.argv[1], int(sys.argv[2]), sys.argv[3].split(",")
    alpha = float(sys.argv[4]) if len(sys.argv) > 4 else 0.35
    method = sys.argv[5] if len(sys.argv) > 5 else "svd"
    order = sys.argv[6] if len(sys.argv) > 6 else "sunfirst"
    stamps, jumps = candidates(hdir, year, bodies, alpha, method, order)
    dates = sorted({s.date().isoformat() for s in stamps})
    print(json.dumps({"year": year, "bodies": [NAMES.get(b, b) for b in bodies], "alpha": alpha,
                      "method": method, "order": order, "candidates": len(stamps),
                      "unique_dates": len(dates), "dates": dates if len(dates) <= 120 else dates[:5] + ["..."],
                      "first_stamps": [s.strftime("%Y-%m-%d %H:%M") for s in stamps[:60]]}, indent=1))


if __name__ == "__main__":
    main()