# Pillar 1: v1.1 (implemented) vs v1.2 (proposed) on synthetic scenarios, 300 seeds each.
import math, random, statistics, sys
GRID = [600 * 72 ** (j / 199) for j in range(200)]
TWO_PI = 2 * math.pi

def qcont(xs, q):
    xs = sorted(xs); h = (len(xs) - 1) * q; lo = math.floor(h); hi = math.ceil(h)
    return xs[lo] + (h - lo) * (xs[hi] - xs[lo])

def core(ts, grid):
    t0 = ts[0]; n = len(ts); Ps = []
    for p in grid:
        f = 1 / p; c = s = 0.0
        for t in ts:
            ph = TWO_PI * f * (t - t0); c += math.cos(ph); s += math.sin(ph)
        Ps.append((c * c + s * s) / n)
    return Ps

def disp(ts):
    d = [b - a for a, b in zip(ts, ts[1:])]
    med = qcont(d, .5)
    if med <= 1e-9: return None, med
    madm = qcont([abs(x - med) for x in d], .5)
    return madm / med, med

def v11(ts):
    ts = sorted(ts)
    if len(ts) < 20: return 0.0, 'insufficient_data', None
    Ps = core(ts, GRID); pmax = max(Ps)
    sig = min(1, max(0, (pmax - math.log(200)) / (6 * math.log(10))))
    dn, med = disp(ts)
    if dn is None: return 0.5 * sig, 'degenerate_deltas', None
    reg = min(1, max(0, 1 - 2 * dn))
    per = max(p for p, v in zip(GRID, Ps) if v >= .5 * pmax)
    return 0.5 * sig + 0.5 * reg, 'ok', per

def v12(ts, cycles=3, gate=0.25):
    ts = sorted(ts)
    if len(ts) < 20: return 0.0, 'insufficient_data', None
    span = ts[-1] - ts[0]
    grid = [p for p in GRID if p <= span / cycles]
    if not grid: return 0.0, 'insufficient_span', None
    Ps = core(ts, grid); pmax = max(Ps); M = len(grid)
    sig = min(1, max(0, (pmax - math.log(M)) / (6 * math.log(10))))
    dn, med = disp(ts)
    if dn is None: return 0.0, 'degenerate_deltas', None
    reg = min(1, max(0, 1 - 2 * dn))
    g = min(1.0, reg / gate)
    per = max(p for p, v in zip(grid, Ps) if v >= .5 * pmax)
    return 0.5 * sig * g + 0.5 * reg, 'ok', per

B = 1_790_000_000.0; T0 = 14400
def perfect(r): return [B + k * T0 for k in range(42)]
def grid(r): return [B + k * T0 + T0 * r.uniform(-.18, .18) for k in range(42)]
def walk(r):
    t, o = B, []
    for _ in range(42): o.append(t); t += T0 * (1 + r.uniform(-.18, .18))
    return o
def poisson(r):
    t, o = B, []
    while len(o) < 40: o.append(t); t += r.expovariate(40 / (7 * 86400))
    return o
def dropout(r):  # 4 h beacon, +-10% cumulative jitter, 30% of beats missed (host asleep)
    t, o = B, []
    for _ in range(60):
        if r.random() > .3: o.append(t)
        t += T0 * (1 + r.uniform(-.1, .1))
    return o
def bursty(r, span):  # browsing-like: Poisson cluster starts, geometric cluster size, ~2 s intra-gaps
    o = []; t = B
    while True:
        t += r.expovariate(1 / 900)
        if t > B + span: break
        k = 1
        while r.random() < .8: k += 1
        u = t
        for _ in range(k): o.append(u); u += r.expovariate(1 / 2)
    return o
def bursty3h(r): return bursty(r, 3 * 3600)
def bursty24h(r): return bursty(r, 86400)
def fast_short(r):  # 15 min beacon observed only 24 h, 10% cumulative jitter
    t, o = B, []
    while t < B + 86400: o.append(t); t += 900 * (1 + r.uniform(-.1, .1))
    return o

def q(xs):
    xs = sorted(xs); n = len(xs)
    return f"min={xs[0]:.3f} p05={xs[n//20]:.3f} med={statistics.median(xs):.3f} p95={xs[-max(1,n//20)]:.3f} max={xs[-1]:.3f}"
seeds = range(1, int(sys.argv[1]) + 1 if len(sys.argv) > 1 else 301)
for name, gen in [('perfect', perfect), ('grid', grid), ('walk', walk), ('dropout', dropout), ('fast_short', fast_short),
                  ('poisson', poisson), ('bursty3h', bursty3h), ('bursty24h', bursty24h)]:
    a, b, pa, pb, st = [], [], [], [], {}
    for s in (seeds if name != 'perfect' else [1]):
        ts = gen(random.Random(s))
        x = v11(ts); y = v12(ts)
        a.append(x[0]); b.append(y[0]); st[y[1]] = st.get(y[1], 0) + 1
        if x[2]: pa.append(x[2])
        if y[2]: pb.append(y[2])
    print(f"{name:10s} v1.1 {q(a)} | >=0.40: {sum(v >= .4 for v in a)}/{len(a)}")
    print(f"{name:10s} v1.2 {q(b)} | >=0.40: {sum(v >= .4 for v in b)}/{len(b)} status={st}")
    if name in ('grid', 'walk', 'dropout'):
        print(f"{'':10s} period within 10% of {T0}: v1.1 {sum(abs(p-T0)<=.1*T0 for p in pa)}/{len(pa)} v1.2 {sum(abs(p-T0)<=.1*T0 for p in pb)}/{len(pb)}")
    if name == 'fast_short':
        print(f"{'':10s} period within 10% of 900: v1.2 {sum(abs(p-900)<=90 for p in pb)}/{len(pb)}")