# C2/C3 floors, independently: min |H| over nonzero frequencies (FFT grid then zoomed grids), the still floor
# by brute minimisation over f' (no closed form), box-k floor by the integer-window argument done numerically
# (min over g real of (g-w)^2 + mu^2 {g}(1-{g}), w = k v), decay floor with mu = mu_s (1-a)/(1+a).
# run: python3 floors_v.py
import numpy as np, vtaps
def hmin(taps, M=1024):
    K = np.zeros((M, M))
    for (dx, dy), w in taps.items(): K[dy % M, dx % M] += w
    A = np.abs(np.fft.fft2(K)); A[0, 0] = np.inf
    # exclude a small neighbourhood of 0 is unnecessary (|H|->1 there)
    j, i = np.unravel_index(np.argmin(A), A.shape)
    fy, fx = ((j + M/2) % M - M/2)/M, ((i + M/2) % M - M/2)/M
    span = 2/M
    for _ in range(12):
        g = np.linspace(-span, span, 41)
        X, Y = np.meshgrid(fx + g, fy + g); V = np.abs(vtaps.H(taps, X, Y))
        k = np.unravel_index(np.argmin(V), V.shape); fx, fy = X[k], Y[k]; span /= 8
    return float(np.abs(vtaps.H(taps, fx, fy))), (fx, fy)
def floor_box(mu, k, nv=4000):
    v = (np.arange(nv) + .5)/nv; w = k*v
    g = np.linspace(-1, k + 1, 20001*k)          # candidate window means x k
    fr = g - np.floor(g)
    best = np.array([np.min((g - wi)**2 + mu*mu*fr*(1 - fr)) for wi in w])/k/k
    return np.sqrt(best.mean())
if __name__ == '__main__':
    for tie in ('up', 'down'):
        t = vtaps.kuw(5, tie); mu, at = hmin(t)
        print(f'tie {tie}: mu_s = {mu:.5f} at {at[0]:.4f},{at[1]:.4f}; still floor {floor_box(mu, 1):.5f}')
    t = vtaps.kuw(5); mu, _ = hmin(t)
    print('box floors:', ' '.join(f'k={k} {floor_box(mu, k):.5f} (still/k {floor_box(mu,1)/k:.5f})' for k in (2, 3, 4, 6)))
    for ms in (20, 40, 80):
        a = np.exp(-1000/60/ms); m = mu*(1 - a)/(1 + a)
        # check min over tau of |G| numerically too
        tau = np.linspace(0, 1, 100001); G = np.abs((1 - a)/(1 - a*np.exp(-2j*np.pi*tau)))
        print(f'decay {ms}: a={a:.5f} min|G|={G.min():.5f} vs (1-a)/(1+a)={(1-a)/(1+a):.5f}; floor {floor_box(m, 1):.5f}')