Commit 666b41ba authored by Philipp Arras's avatar Philipp Arras
Browse files

Start working on proper w-stacking tests

parent 1e121c2b
......@@ -16,7 +16,6 @@ def _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow):
speedoflight, f0 = 3e8, 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight)
uvw[:, 2] = 0.
baselines = ng.Baselines(coord=uvw, freq=freq)
flags = np.zeros((nrow, nchan), dtype=np.bool)
idx = ng.getIndices(baselines, conf, flags)
......@@ -260,3 +259,82 @@ def test_against_dft(nxdirty, nydirty, epsilon, nchan, nrow):
phase = x*uvw[ii, 0] + y*uvw[ii, 1]
res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
assert_(_l2error(res0, res1) < epsilon)
@pmp('epsilon', [1e-7])
@pmp('nxdirty', [12, 64])
@pmp('nydirty', [64])
@pmp("nrow", (10, 100, 1000, 10000))
@pmp("nchan", (1,))
def test_against_wdft(nxdirty, nydirty, epsilon, nchan, nrow):
np.random.seed(40)
pixsize = np.pi/180/60/nxdirty
conf = ng.GridderConfig(nxdirty=nxdirty,
nydirty=nydirty,
epsilon=epsilon,
pixsize_x=pixsize,
pixsize_y=pixsize)
speedoflight, f0 = 3e8, 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight)
bl = ng.Baselines(coord=uvw, freq=freq)
flags = np.zeros((nrow, nchan), dtype=np.bool)
idx = ng.getIndices(bl, conf, flags)
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
vis = bl.ms2vis(ms, idx)
uvw = bl.effectiveuvw(idx)
res0 = np.zeros((nxdirty, nydirty))
mi, ma = np.min(uvw[:, 2]), np.max(uvw[:, 2])
ws = np.linspace(mi, ma, 2)
delta = (ma-mi)/10000
ws = np.arange(mi, ma, delta)
for ww in ws:
jj = ng.getIndices(bl, conf, flags, wmin=ww, wmax=ww+delta)
if len(jj) == 0:
continue
dd = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, jj, bl.ms2vis(ms, jj)))
res0 += conf.apply_wscreen(dd, ww, adjoint=True).real
# Compute dft with w term
x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty]],
indexing='ij')
x *= conf.Pixsize_x()
y *= conf.Pixsize_y()
res1 = np.zeros_like(res0)
n = np.cos(np.sqrt(x**2 + y**2))
for ii in idx:
phase = x*uvw[ii, 0] + y*uvw[ii, 1] + uvw[ii, 2]*(n-1)
print(np.max(x*uvw[ii, 0]), np.max(y*uvw[ii, 1]),np.max(uvw[ii, 2]*(n-1)) )
res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
res1 /= n
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(2, 2, 1)
ax.set_title('wstacking')
im = ax.imshow(res0)
fig.colorbar(im, ax=ax)
ax = fig.add_subplot(2, 2, 2)
ax.set_title('ground truth')
im = ax.imshow(res1)
fig.colorbar(im, ax=ax)
ax = fig.add_subplot(2, 2, 3)
ax.set_title('relative error')
im = ax.imshow(np.abs(res1-res0)/res1)
fig.colorbar(im, ax=ax)
ax = fig.add_subplot(2, 2, 4)
ax.set_title('error, l2 = {}'.format(_l2error(res0, res1)))
im = ax.imshow(res1-res0)
fig.colorbar(im, ax=ax)
plt.savefig('nx{}_ny{}_nchan{}_nrow{}_eps{}.png'.format(nxdirty, nydirty, nchan, nrow, epsilon))
plt.close()
assert_(_l2error(res0, res1) < epsilon)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment