# Paint off, one pixel alone (the identity in space): how well can a 0/1 sequence over frames, blended by the
# eye, hold a level f? A box over k frames: exactly dist(f, multiples of 1/k), RMS over levels 1/(sqrt(12) k).
# A decay y_t = a y_{t-1} + (1 - a) b_t: an average-cost control problem, solved by relative value iteration on
# a grid of y; the policy it finds is then run, and what it achieves is reported too (so the number is
# attained, not only bounded). Also the Parseval floor ((1 - a)/(1 + a) is the least a decay passes), a greedy
# "leaky sigma-delta" (pick the bit that lands nearer f), and the best one-pixel lattice C/N.
# (From the theory pass of this run, with the eye's time constants matched to st.py.)
import numpy as np

def parseval_floor(f, g):
    fp = np.clip((2 * f - g * g) / (2 - 2 * g * g), 0, 1)
    return g * g * fp * (1 - fp) + (f - fp) ** 2

def dp_opt(f, a, M=4001, iters=2500, sim=False):
    """The least long-run mean of (y_t - f)^2 over 0/1 sequences, by relative value iteration (y on M points,
    linear interpolation). With sim, also the mean the greedy-on-V policy achieves over 40000 frames."""
    y = np.linspace(0, 1, M)
    def interp(V, yq):
        p = yq * (M - 1); i = np.clip(np.floor(p).astype(int), 0, M - 2); w = p - i
        return V[i] * (1 - w) + V[i + 1] * w
    y0 = a * y; y1 = a * y + (1 - a)
    c0 = (y0 - f) ** 2; c1 = (y1 - f) ** 2
    V = np.zeros(M); lams = []; ref = int(round(f * (M - 1)))
    for _ in range(iters):
        Vn = np.minimum(c0 + interp(V, y0), c1 + interp(V, y1))
        lams.append(Vn[ref]); V = Vn - Vn[ref]
    lam = float(np.mean(lams[-iters // 4:]))   # (a deterministic problem can cycle: average the increments)
    if not sim: return lam
    pol = c1 + interp(V, y1) < c0 + interp(V, y0)
    yy = f; e2 = 0.0; T = 40000
    for t in range(T):
        yy = a * yy + (1 - a) * (1.0 if pol[int(round(yy * (M - 1)))] else 0.0)
        if t >= 2000: e2 += (yy - f) ** 2
    return lam, e2 / (T - 2000)

def greedy(f, a, T=60000):
    y = f; e2 = 0.0
    for t in range(T):
        y0 = a * y; y1 = y0 + 1 - a
        y = y1 if abs(y1 - f) < abs(y0 - f) else y0
        if t >= 1000: e2 += (y - f) ** 2
    return e2 / (T - 1000)

def lattice1d(a, N, C):
    m = np.arange(1, N); tau = m * C / N
    G2 = (1 - a) ** 2 / (1 - 2 * a * np.cos(2 * np.pi * tau) + a * a)
    return 1 / (12 * N * N) + (G2 / (2 * N * N * np.sin(np.pi * m / N) ** 2)).sum()

if __name__ == '__main__':
    print('paint off, one pixel: RMS over 32 levels')
    print('  box k (exact): ' + ', '.join(f'k={k} {1 / np.sqrt(12) / k:.4f}' for k in (2, 3, 4, 6)))
    fav = (np.arange(32) + .5) / 32
    for ms in (20, 40, 80):
        a = np.exp(-1000 / 60 / ms); g = (1 - a) / (1 + a)
        P = np.mean([parseval_floor(f, g) for f in fav]); DD = [dp_opt(f, a, sim=True) for f in fav]
        D = np.mean([d for d, _ in DD]); Ds = np.mean([s for _, s in DD]); Gr = np.mean([greedy(f, a) for f in fav])
        best = min((lattice1d(a, N, C), C, N) for N in range(2, 65) for C in range(1, N) if np.gcd(C, N) == 1)
        print(f'  decay {ms} ms: Parseval floor {np.sqrt(P):.4f}, the best sequence {np.sqrt(D):.4f} (its policy run: {np.sqrt(Ds):.4f}), '
              f'greedy {np.sqrt(Gr):.4f}, best lattice C/N = {best[1]}/{best[2]} {np.sqrt(best[0]):.4f}, flip every frame {np.sqrt(lattice1d(a, 2, 1)):.4f}, white {np.sqrt(g / 6):.4f}')