#!/usr/bin/env python3
# B7A analysis: offsets of DMS / NITZ / NTP(vniiftri 1..3) vs the board's free-running CLOCK_REALTIME (1970-based), then DMS/NITZ - NTP.
import re, sys, datetime, statistics as st
L = open(sys.argv[1]).read().splitlines(); GPS = 315964800
dms, nitz, ntp = [], [], []
for l in L:
    m = re.match(r'r(\d) DMS\d up=([\d.]+) t0=([\d.]+) t1=([\d.]+) (\d+)$', l)
    if m: r,u,a,b,ms = int(m[1]),float(m[2]),float(m[3]),float(m[4]),int(m[5]); dms.append((r,u,(a+b)/2,(b-a)/2, ms/1000+GPS-(a+b)/2)); continue
    m = re.match(r'r(\d) NITZ\d up=([\d.]+) t0=([\d.]+) t1=([\d.]+) (\S+) (\S+) tz=(\d+)', l)
    if m:
        r,u,a,b = int(m[1]),float(m[2]),float(m[3]),float(m[4])
        e = datetime.datetime.fromisoformat(m[5]+'T'+m[6]+'+00:00').timestamp(); nitz.append((r,u,(a+b)/2,(b-a)/2,e-(a+b)/2)); continue
    m = re.match(r'r(\d) NTP(\d) up=([\d.]+) t0=([\d.]+) t1=([\d.]+)', l)
    if m: cur = (int(m[1]),int(m[2]),float(m[3]),float(m[4])); continue
    m = re.match(r'r(\d) NTP(\d) ntpd: reply from \S+: offset:([+-][\d.]+) delay:([\d.]+)', l)
    if m: ntp.append((int(m[1]),int(m[2]),float(m[3]),float(m[4])))
# NTP: best (min delay) reply per round/server; local time of the run start t0 approximates sample time (8 s window)
starts = {}
for l in L:
    m = re.match(r'r(\d) NTP(\d) up=([\d.]+) t0=([\d.]+) t1=([\d.]+)', l)
    if m: starts[(int(m[1]),int(m[2]))] = (float(m[4])+float(m[5]))/2
best = {}
for r,s,o,d in ntp:
    if (r,s) not in best or d < best[(r,s)][1]: best[(r,s)] = (o,d)
# linear fit of NTP offset vs local time (all best samples, all servers) -> local clock drift
xs = [starts[k] for k in best]; ys = [best[k][0] for k in best]
mx, my = st.mean(xs), st.mean(ys); slope = sum((x-mx)*(y-my) for x,y in zip(xs,ys))/sum((x-mx)**2 for x in xs)
ntp_at = lambda t: my + slope*(t-mx)
print(f"NTP samples: {len(ntp)} replies, best per round/server: {len(best)}; local clock drift vs NTP: {slope*1e6:+.2f} ppm")
print("\nPer-round NTP best offsets (s, relative to board clock) and server spread:")
print("round  ntp1(delay)                 ntp2(delay)                 ntp3(delay)                 spread_ms  resid_vs_fit_ms")
for r in sorted({k[0] for k in best}):
    row = [best.get((r,s)) for s in (1,2,3)]
    os_ = [x[0] for x in row if x]
    res = [ (best[(r,s)][0]-ntp_at(starts[(r,s)]))*1e3 for s in (1,2,3) if (r,s) in best]
    print(f"r{r}  " + "  ".join(f"{x[0]:.6f}({x[1]*1e3:5.1f}ms)" if x else "-" for x in row) + f"  {(max(os_)-min(os_))*1e3:7.2f}  " + " ".join(f"{v:+.2f}" for v in res))
print("\nTable: uptime / DMS / NITZ / ntp1 / ntp2 / ntp3 (UTC at the sample's local midpoint; NTP = board clock + NTP offset at that instant)")
def utc(t): return datetime.datetime.fromtimestamp(t, datetime.timezone.utc).strftime('%H:%M:%S.%f')[:-3]
print("round uptime   DMS_UTC        NITZ_UTC   ntp1           ntp2           ntp3           DMS-NTP_ms  qmi_halfwin_ms")
dd = []
for r in sorted({x[0] for x in dms}):
    d = [x for x in dms if x[0]==r]; n = [x for x in nitz if x[0]==r]
    for x in d:
        t = x[2]; diff = (x[4]-ntp_at(t))*1e3; dd.append((diff, x[3]*1e3))
        no = [best[(r,s)][0] for s in (1,2,3) if (r,s) in best]
        cols = [utc(t+best[(r,s)][0]) if (r,s) in best else '-' for s in (1,2,3)]
        print(f"r{r}   {x[1]:8.2f}  {utc(t+x[4])}  {utc(t+n[0][4])[:8] if n else '-':10} {cols[0]}  {cols[1]}  {cols[2]}  {diff:+8.1f}    {x[3]*1e3:5.1f}")
nd = [((x[4]-ntp_at(x[2]))*1e3, x[3]*1e3) for x in nitz]
print("\nDMS - NTP (ms): n=%d mean %+.1f median %+.1f min %+.1f max %+.1f stdev %.1f ; max |DMS-NTP| %.1f ms; qmicli half-window max %.1f ms" % (
    len(dd), st.mean(v for v,_ in dd), st.median(v for v,_ in dd), min(v for v,_ in dd), max(v for v,_ in dd), st.pstdev([v for v,_ in dd]), max(abs(v) for v,_ in dd), max(h for _,h in dd)))
print("NITZ - NTP (ms, NITZ truncated to 1 s => expected in (-1000, 0]): n=%d min %+.1f max %+.1f" % (len(nd), min(v for v,_ in nd), max(v for v,_ in nd)))
for x in nitz: pass
print("NTP server spread (max-min of best offsets within a round): max %.2f ms" % max((max(best[(r,s)][0] for s in (1,2,3) if (r,s) in best)-min(best[(r,s)][0] for s in (1,2,3) if (r,s) in best))*1e3 for r in {k[0] for k in best}))
print("DMS formula check: Unix UTC = DMS_ms/1000 + 315964800 -> DMS-NTP mean above (0 = formula exact; ~ +18000 ms would mean GPS time without leap correction)")