SSM-JTK – CSV-only Scorecard How-To (5)

# scorecard_csv_only.py (plain ASCII; reads two CSVs with identical date/planet coverage)
# Usage: python scorecard_csv_only.py golden.csv reference.csv Venus
import sys, csv, math, datetime as dt

def wrap360(x): return (x % 360.0 + 360.0) % 360.0
def wrap180(d):
    y = (d + 180.0) % 360.0
    return y - 180.0
def rasi(x): return int(math.floor(wrap360(x) / 30.0))
def cusp_dist(x):
    q = wrap360(x) % 30.0
    return min(q, 30.0 - q)

def read_csv(path):
    rows = []
    with open(path, newline="", encoding="utf-8") as f:
        for r in csv.DictReader(f):
            rows.append({
                "planet": r["planet"].strip().lower(),
                "date": dt.date.fromisoformat(r["date"]),
                "L": float(r.get("lon_sidereal_lahiri_deg") or r.get("L_actual_deg"))
            })
    return rows

def group_by(rows, planet):
    return [r for r in rows if r["planet"] == planet.lower()]

def unwrap(series):
    y = [series[0]]
    for i in range(1, len(series)):
        x = series[i]; best = x
        for k in (-1, 0, 1):
            cand = x + 360.0 * k
            if abs(cand - y[-1]) < abs(best - y[-1]): best = cand
        y.append(best)
    return y

def central_diff(y):
    if len(y) < 3: return [0.0]*len(y)
    return [y[1]-y[0]] + [(y[i+1]-y[i-1])/2.0 for i in range(1,len(y)-1)] + [y[-1]-y[-2]]

def detect_crossings(t, y_unwrapped, step=30.0):
    cr = []
    y0, y1 = min(y_unwrapped[0], y_unwrapped[-1]), max(y_unwrapped[0], y_unwrapped[-1])
    kmin = math.floor(y0/step) - 1; kmax = math.ceil(y1/step) + 1
    for k in range(int(kmin), int(kmax)+1):
        target = k*step
        for i in range(1, len(y_unwrapped)):
            s0 = y_unwrapped[i-1] - target
            s1 = y_unwrapped[i]   - target
            if s0 == 0.0: cr.append(t[i-1]); continue
            if s0*s1 < 0.0:
                frac = (target - y_unwrapped[i-1]) / (y_unwrapped[i] - y_unwrapped[i-1])
                tau  = t[i-1] + frac*(t[i]-t[i-1])
                cr.append(tau)
    return cr

def pair_events(A, B, cap_days):
    used, pairs = set(), []
    for a in A:
        jbest, dbest = None, None
        for j,b in enumerate(B):
            if j in used: continue
            d = abs(b - a)
            if dbest is None or d < dbest:
                jbest, dbest = j, d
        if dbest is not None and dbest <= cap_days:
            pairs.append((a, B[jbest])); used.add(jbest)
    return pairs

def main(golden_path, reference_path, planet):
    G = group_by(read_csv(golden_path), planet)
    R = group_by(read_csv(reference_path), planet)
    if len(G) != len(R):
        raise SystemExit("Mismatched coverage; align date/planet first.")

    G.sort(key=lambda r: r["date"]); R.sort(key=lambda r: r["date"])
    dates = [r["date"] for r in G]
    Lg = [r["L"] for r in G]; Lr = [r["L"] for r in R]

    # pointwise degree metrics
    errs = [abs(wrap180(Lg[i]-Lr[i])) for i in range(len(dates))]
    mae = sum(errs)/len(errs)
    p90 = sorted(errs)[int(0.90*len(errs))-1]
    mxx = max(errs)
    mis = sum(1 for i in range(len(dates)) if rasi(Lg[i])!=rasi(Lr[i]))/len(dates)

    # cusp MAE
    cusp_mae = sum(abs(cusp_dist(Lg[i]) - cusp_dist(Lr[i])) for i in range(len(dates))) / len(dates)

    # events (unwrapped)
    t   = [(d - dates[0]).days for d in dates]
    yG  = unwrap(Lg); yR = unwrap(Lr)

    # crossings
    cG  = detect_crossings(t, yG, 30.0)
    cR  = detect_crossings(t, yR, 30.0)
    Cp  = pair_events(cR, cG, cap_days=60.0)
    cross_mae = (sum(abs(a-b) for a,b in Cp)/len(Cp)) if Cp else float("nan")

    # stations (simple)
    vG, vR = central_diff(yG), central_diff(yR)
    sG = [t[i] for i in range(1,len(vG)-1) if abs(vG[i])<=abs(vG[i-1]) and abs(vG[i])<=abs(vG[i+1])]
    sR = [t[i] for i in range(1,len(vR)-1) if abs(vR[i])<=abs(vR[i-1]) and abs(vR[i])<=abs(vR[i+1])]
    Sp = pair_events(sR, sG, cap_days=90.0)
    stat_mae = (sum(abs(a-b) for a,b in Sp)/len(Sp)) if Sp else float("nan")

    print(f"Planet={planet}")
    print(f"MAE_deg={mae:.6f}  P90_deg={p90:.6f}  MAX_deg={mxx:.6f}  misclass_rate={mis:.6f}")
    print(f"cusp_dist_MAE_deg={cusp_mae:.6f}")
    print(f"rasi_cross_MAE_days={cross_mae:.6f}")
    print(f"station_date_MAE_days={stat_mae:.6f}")

if __name__ == "__main__":
    if len(sys.argv) != 4:
        print("Usage: python scorecard_csv_only.py golden.csv reference.csv PlanetName")
        sys.exit(2)
    main(sys.argv[1], sys.argv[2], sys.argv[3])