kernel_helper.py 3.77 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
import numpy as np

def gridder_to_C(gridder, W):
    M = len(gridder) // W
    C = np.zeros((W, M), dtype=float)
    for r in range(0, W):
        ell = r - (W / 2) + 1
        indx = (np.arange(M) - 2 * M * ell).astype(int)
        # Use symmetry to deal with negative indices
        indx[indx < 0] = -indx[indx < 0] - 1
        C[r, :] = gridder[indx]
    return C

def C_to_grid_correction(C, nu_C, x, optimal=True):
    W = C.shape[0]
    nx = x.shape[0]
    c = np.zeros(x.shape, dtype=float)
    d = np.zeros(x.shape, dtype=float)
    C = C.reshape((C.shape[0],C.shape[1],1))
    nu_C = nu_C.reshape((-1,1))
    x = x.reshape((1,-1))
    cosarr = np.empty((W,nx))
    for i in range(W):
        cosarr[i,:] = np.cos(2 * np.pi * i * x)
    for rp in range(0, W):
        ellp = rp - (W / 2) + 1
        for r in range(0, W):
            ell = r - (W / 2) + 1
            xmn = np.mean(C[rp, :, :]*C[r, :, :], axis=0)
            tmp = xmn * cosarr[abs(rp-r)]
            d += tmp
        tmp2 = C[rp, :, :] * np.cos(2 * np.pi * (ellp-nu_C) * x)
        c += np.mean(tmp2,axis=0)
    return c / d if optimal else 1 / c


def gridder_to_grid_correction(gridder, nu, x, W, optimal=True):
    M = len(nu) // W
    C = gridder_to_C(gridder, W)
    return C_to_grid_correction(C, nu[:M], x, optimal)


def calc_map_error_from_C(C, grid_correction, nu_C, x, W):
    M = len(nu_C)
    nx = len(x)
    one_app=np.zeros((nx, M), dtype=np.complex128)
    for r in range(0, W):
        ell = r - (W / 2) + 1
        one_app += grid_correction.reshape((-1,1)) * C[r, :] \
            * np.exp(2j * np.pi * (ell - nu_C).reshape((1,-1)) * x.reshape((-1,1)))
    one_app = (1.-one_app.real)**2 + one_app.imag**2
    map_error = np.sum(one_app,  axis=1)/M
    return map_error


def calc_map_error(gridder, grid_correction, nu, x, W):
    M = len(nu) // W
    C = gridder_to_C(gridder, W)
    return calc_map_error_from_C(C, grid_correction, nu[:M], x, W)


def eskapprox(parm, nu, x, W):
    nunorm=2*nu/W
    beta=parm[0]
    e1 = 0.5 if len(parm)<2 else parm[1]
    e2 = 2. if len(parm)<3 else parm[2]
    return np.exp(beta*W*((1-nunorm**e2)**e1-1))

def getmaxerr(approx, coeff, nu, x, W, M, N, x0):
    nu=(np.arange(W*M)+0.5)/(2*M)
    x=np.arange(N+1)/(2*N)
    krn = approx(coeff, nu, x, W)
    err = kernel2error(krn, nu, x, W)
    err = err[0:int(2*x0*N+0.9999)+1]
    return np.max(np.abs(err))

Martin Reinecke's avatar
Martin Reinecke committed
77
def scan_esk(rbeta, re0, nu, x, W, M, N, x0, nsamp):
Martin Reinecke's avatar
Martin Reinecke committed
78
    curmin=1e30
Martin Reinecke's avatar
Martin Reinecke committed
79
80
    for e0 in np.linspace(re0[0], re0[1], nsamp):
        for beta in np.linspace(rbeta[0], rbeta[1], nsamp):
Martin Reinecke's avatar
Martin Reinecke committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
            test = getmaxerr(eskapprox, [beta,e0], nu, x, W, M, N, x0)
            if test<curmin:
                curmin, coeffmin = test, [beta,e0]
    return coeffmin

def kernel2error(krn, nu, x, W):
    corr = gridder_to_grid_correction(krn, nu, x, W)
    return calc_map_error(krn, corr, nu, x, W)

M=128
N=512
x=np.arange(N+1)/(2*N)

# for quick experiments, just enter the desired oversampling factor and support
# as single elements in the tuples below
Martin Reinecke's avatar
Martin Reinecke committed
96
97
ofactors = np.linspace(1.15,2.00,18)
Ws = np.arange(4,17)
Martin Reinecke's avatar
Martin Reinecke committed
98
99
100
101
102
103
104
results = []
for W in Ws:
    for ofactor in ofactors:
        x0 = 0.5/ofactor
        nu=(np.arange(W*M)+0.5)/(2*M)
        ulim = int(2*x0*N+0.9999)+1
        rbeta=[1., 2.5]
Martin Reinecke's avatar
Martin Reinecke committed
105
        re0=[0.48, 0.65]
Martin Reinecke's avatar
Martin Reinecke committed
106
107
        dbeta = rbeta[1]-rbeta[0]
        de0 = re0[1]-re0[0]
Martin Reinecke's avatar
Martin Reinecke committed
108
109
110
111
        for i in range(30):
            res1 = scan_esk(rbeta, re0, nu, x, W, M, N, x0, 10)
            dbeta*=0.5
            de0*=0.5
Martin Reinecke's avatar
Martin Reinecke committed
112
113
114
115
116
            rbeta = [res1[0]-0.5*dbeta, res1[0]+0.5*dbeta]
            re0 = [res1[1]-0.5*de0, res1[1]+0.5*de0]
        krn1 = eskapprox(res1, nu, x, W) 
        err1 = kernel2error(krn1, nu, x, W)
        maxerr1 = np.sqrt(np.max(err1[0:ulim]))
Martin Reinecke's avatar
Martin Reinecke committed
117
        print("{{{0:2d}, {1:4.2f}, {2:13.8g}, {3:12.10f}, {4:12.10f}}},".format(W, ofactor, maxerr1, res1[0], res1[1]))