# C4 (golden shift) and C5 (practical dithers) by my own simulation: round at 32 flat levels, blend over frames
# FIRST (causal box / decay with a long warm-up), then filter in space by direct shifted sums on a canvas with a
# margin (no FFT, no wrap). White-noise theory from sum h^2 and sum g^2.  run: python3 c5_sim.py [x0 y0]
import sys, math, numpy as np, vtaps
taps = vtaps.kuw(5); R = 5; KEEP = 96; n = KEEP + 2*R
X0, Y0 = (float(sys.argv[1]), float(sys.argv[2])) if len(sys.argv) > 2 else (60., 1100.)
PHI = (math.sqrt(5) - 1)/2
WARM, TK = 160, 64; T = WARM + TK
f32 = np.float32
def fr32(v): return (v - np.floor(v)).astype(f32)
def ign(x, y):
    x = x.astype(f32); y = y.astype(f32)
    return fr32(f32(52.9829189)*fr32(x*f32(0.06711056) + y*f32(0.00583715))).astype(float)
yy, xx = np.mgrid[0:n, 0:n] + .5; xx = xx + X0; yy = yy + Y0
def frames(name):
    out = np.empty((T, n, n)); rng = np.random.default_rng(3)
    base = ign(xx, yy)
    lat = (((54*(xx - .5).astype(int) + 11*(yy - .5).astype(int)) % 56) + .5)/56
    for t in range(T):
        if name == 'ign still': out[t] = base
        elif name == 'ign half': out[t] = (base + .5*(t % 2)) % 1
        elif name == 'ign golden': out[t] = (base + PHI*t) % 1
        elif name == 'ign wolfe': s = 5.588238*(t % 64); out[t] = ign(xx + s, yy + s)
        elif name == 'white': out[t] = rng.random((n, n))
        elif name == 'lattice golden': out[t] = (lat + PHI*t) % 1
    return out
def eyes():
    out = []
    for k in (2, 4): out.append((f'box {k}', ('box', k)))
    for ms in (40, 80): out.append((f'decay {ms}', ('decay', math.exp(-1000/60/ms))))
    return out
def blend(e, eye):
    if eye[0] == 'box':
        k = eye[1]; return sum(e[WARM - s:T - s] for s in range(k))/k
    a = eye[1]; y = np.zeros_like(e[0]); out = []
    for t in range(T):
        y = a*y + (1 - a)*e[t]
        if t >= WARM: out.append(y.copy())
    return np.array(out)
def spatial(img):
    out = np.zeros(img.shape[:-2] + (KEEP, KEEP))
    for (dx, dy), w in taps.items(): out += w*img[..., R + dy:R + dy + KEEP, R + dx:R + dx + KEEP]
    return out
def score(name, levels=32):
    fr_ = frames(name); acc = np.zeros(4)
    for v in (np.arange(levels) + .5)/levels:
        e = np.floor(v + fr_) - v
        for i, (_, eye) in enumerate(eyes()): acc[i] += (spatial(blend(e, eye))**2).mean()
    return np.sqrt(acc/levels)
if __name__ == '__main__':
    print(f'origin ({X0},{Y0}); columns: ' + ' '.join(nm for nm, _ in eyes()))
    for nm in ['ign still', 'ign half', 'ign golden', 'ign wolfe', 'white', 'lattice golden']:
        print(f'  {nm:15s} ' + ' '.join(f'{x:.4f}' for x in score(nm)), flush=True)
    h2 = sum(w*w for w in taps.values()); base = math.sqrt(h2/6)
    g2 = [.5, .25] + [(1 - a)/(1 + a) for a in (math.exp(-1000/60/40), math.exp(-1000/60/80))]
    print(f'  white theory    ' + ' '.join(f'{base*math.sqrt(g):.4f}' for g in g2) + f'   (sqrt(sum h^2/6) = {base:.4f})')
    f = frames('ign half'); print('  half-flip period 2 exactly:', bool(np.all(f[2:] == f[:-2])))