# Floors for any dither in space and time (any pattern, any rule that outputs integers), under each eye model,
# against the best lattice found (search.py). Prints out/floors.txt.
#
# A box over k frames: the sum of the k frames a pixel shows is an integer image S in {0 .. k}, and the eye sees
# h * S / k. If S has mean k f', Parseval gives (f' - f)^2 + sum over nonzero frequencies of |H|^2 |S^|^2 / k^2, and
# an integer image with mean k f' has variance at least {k f'}(1 - {k f'}). So at level f the mean-square error is
# at least min over f' of (f' - f)^2 + mu^2 {k f'}(1 - {k f'}) / k^2 (mu the smallest |H| at a nonzero spatial
# frequency): the still floor with the level grid 1/k finer, so the RMS floor is the still one divided by k.
# A decay: the eye's blend is a linear filter in space and time with response H(xi) G(tau), and the frames are
# integer images, so the still floor's argument runs on the space-time image with mu = min |H G| over nonzero
# (xi, tau): mu_s (1 - a) / (1 + a), at tau = 1/2 (a cosine flipping every frame is what a decay passes least).
# A decay, stronger (an estimate, not certified): per frame, Parseval in space gives the error at least
# (f'_t - f)^2 + mu_s^2 Var_p(y_t), and each pixel's blended value y_t over time can hold its mean no better than
# the best 0/1 sequence does (temporal1d.py, numerically: E1(f')); so E(f) >= min over f' of (f' - f)^2 + mu_s^2 E1(f').
import numpy as np
import filters as F, lattice as L, st
from temporal1d import dp_opt

def box_floor(still_rms, k):
    return still_rms / k

taps = F.kuwahara(5)
bd = L.bounds(taps)
mu = bd['hmin']
print(f'Starlock full pass: smallest |H| at a nonzero spatial frequency {mu:.4f}; still floor {bd["any"]:.4f}')
# (the best lattice per eye model, from out/search.txt)
BEST = {'box 2': (25, 5, 26, 52), 'box 3': (20, 4, 21, 63), 'box 4': (1, 13, 16, 64), 'box 6': (59, 20, 10, 61),
        'decay 20 ms': (55, 11, 34, 57), 'decay 40 ms': (10, 2, 21, 51), 'decay 80 ms': (12, 2, 39, 63)}
print(f'  {"eye":12s} {"floor":>7s} {"best lattice":>18s} {"RMS":>7s} {"x floor":>8s}')
for G in [st.box(2), st.box(3), st.box(4), st.box(6), st.decay(20), st.decay(40), st.decay(80)]:
    if G.taps is not None: fl = box_floor(bd['any'], len(G.taps)); est = ''
    else:
        fl = L.any_floor(mu * (1 - G.a) / (1 + G.a))
        fp = np.linspace(0, 1, 201); E1 = np.array([dp_opt(x, G.a) for x in fp]); f = (np.arange(4000) + .5) / 4000
        e = np.sqrt(((fp[None] - f[:, None]) ** 2 + mu * mu * E1[None]).min(1).mean())
    A, B, C, N = BEST[G.name]; j = np.sqrt(st.J(A, B, C, N, taps, G))
    print(f'  {G.name:12s} {fl:7.4f} {f"({A},{B},{C})/{N}":>18s} {j:7.4f} {j / fl:8.3f}' + ('' if G.taps is not None else f'   (stronger estimate {e:.4f}: {j / e:.3f} x)'))