Skip to content

GitLab

Commit 000144b1 authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'sometests' into separate_interfaces

parents 3e52d27a 72de39d7
Hide whitespace changes
Inline Side-by-side
test.py
View file @ 000144b1
......@@ -11,8 +11,8 @@ def _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow):
conf = ng.GridderConfig(nxdirty=nxdirty,
nydirty=nydirty,
epsilon=epsilon,
pixsize_x=pixsize,
pixsize_y=0.5*pixsize)
pixsize_x=0.368*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)
......@@ -22,12 +22,24 @@ def _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow):
return baselines, conf, idx
def _wscreen(npix, dst, w):
dc = (np.linspace(start=-npix/2, stop=npix/2 - 1, num=npix)*dst)**2
ls = np.broadcast_to(dc, (dc.shape[0],)*2)
n = np.sqrt(1-ls-ls.T)
wscreen = np.exp(2*np.pi*1j*w*(n - 1))/n
return wscreen
def _dft(uvw, vis, conf, apply_w):
shp = (conf.Nxdirty(), conf.Nydirty())
x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in shp],
indexing='ij')
x *= conf.Pixsize_x()
y *= conf.Pixsize_y()
dirty = np.zeros(shp, dtype=np.complex128)
if apply_w:
n = np.sqrt(1-x**2-y**2)
nm1 = n-1
for ii, vv in enumerate(vis):
phase = x*uvw[ii, 0] + y*uvw[ii, 1]
if apply_w:
phase += uvw[ii, 2]*nm1
dirty += (vv*np.exp(2j*np.pi*phase))
if apply_w:
return dirty/n
return dirty
def _l2error(a, b):
......@@ -47,10 +59,26 @@ def test_adjointness(nxdirty, nydirty, nrow, nchan, epsilon):
grid = ng.vis2grid_c(bl, gconf, idx, vis)
dirty = gconf.grid2dirty_c(grid)
dirty2 = np.random.rand(*dirty.shape)-0.5
ms2 = bl.vis2ms(ng.grid2vis_c(bl, gconf, idx, gconf.dirty2grid_c(dirty2)), idx)
ms2 = bl.vis2ms(ng.grid2vis_c(bl, gconf, idx, gconf.dirty2grid_c(dirty2)),
idx)
assert_allclose(np.vdot(ms, ms2).real, np.vdot(dirty.real, dirty2.real))
@pmp("nxdirty", (128, 300))
@pmp("nydirty", (128, 250))
@pmp("nrow", (1, 10, 10000))
@pmp("nchan", (1, 10, 100))
@pmp("epsilon", (1e-2, 1e-7, 2e-13))
def test_adjointness_wgridding(nxdirty, nydirty, nrow, nchan, epsilon):
np.random.seed(42)
bl, conf, idx = _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow)
vis = np.random.rand(*idx.shape)-0.5
dirty = np.random.rand(conf.Nxdirty(), conf.Nydirty())-0.5
dirty2 = ng.vis2dirty_wstack(bl, conf, idx, vis)
vis2 = ng.dirty2vis_wstack(bl, conf, idx, dirty)
assert_allclose(np.vdot(vis, vis2), np.vdot(dirty2, dirty), rtol=epsilon)
@pmp("nxdirty", (128,))
@pmp("nydirty", (128,))
@pmp("nrow", (10000,))
......@@ -137,26 +165,6 @@ def test_pickling():
assert_(gc.Pixsize_y() == gc2.Pixsize_y())
@pmp('nx', [4, 18, 54])
@pmp('fov', [0.1, 1, 5]) # deg
@pmp('w', [0, 10, 8489])
def test_wstacking(nx, fov, w):
np.random.seed(42)
dx = fov*np.pi/180/nx
ny, dy = nx, dx
conf = ng.GridderConfig(nx, ny, 1e-7, dx, dy)
x, y = conf.Nu(), conf.Nv()
fld = np.random.randn(x, y) + 1j*np.random.randn(x, y)
res0 = (_wscreen(nx, dy, w).conjugate()*conf.grid2dirty_c(fld)).real
res1 = conf.apply_wscreen(conf.grid2dirty_c(fld), w, True).real
np.testing.assert_allclose(res0, res1)
fld = np.random.randn(nx, ny)
res0 = conf.dirty2grid_c(_wscreen(nx, dx, w)*fld.real)
res1 = conf.dirty2grid_c(conf.apply_wscreen(fld, w, False))
np.testing.assert_allclose(res0, res1)
@pmp("nxdirty", (128, 300))
@pmp("nydirty", (128, 250))
@pmp("nrow", (1, 10, 10000))
......@@ -205,6 +213,24 @@ def test_correlations(nxdirty, nydirty, nrow, nchan, epsilon, du, dv, weight):
assert_allclose(np.zeros_like(y1), y1.imag)
@pmp("nrow", (1, 1000))
@pmp("nchan", (1, 10))
@pmp("epsilon", (1e-2, 1e-7, 2e-12))
@pmp("weight", (True, False))
def test_no_correlation(nrow, nchan, epsilon, weight):
np.random.seed(42)
bl, conf, idx = _init_gridder(128, 128, epsilon, nchan, nrow)
wgt = np.random.rand(*idx.shape) if weight else None
w = conf.Supp()
for uu in range(-2*w, 2*w):
for vv in range(-2*w, 2*w):
if abs(uu) >= w or abs(vv) >= w:
with pytest.raises(RuntimeError):
ng.get_correlations(bl, conf, idx, uu, vv, wgt)
else:
ng.get_correlations(bl, conf, idx, uu, vv, wgt)
@pmp('epsilon', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13])
@pmp('nxdirty', [12, 128])
@pmp('nydirty', [4, 12, 128])
......@@ -214,18 +240,9 @@ def test_against_dft(nxdirty, nydirty, epsilon, nchan, nrow):
bl, conf, idx = _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow)
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
vis = bl.ms2vis(ms, idx)
res0 = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, idx, vis))
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)
uvw = bl.effectiveuvw(idx)
for ii in idx:
phase = x*uvw[ii, 0] + y*uvw[ii, 1]
res1 += (vis[ii]*np.exp(2j*np.pi*phase))
res0 = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, idx, vis))
res1 = _dft(uvw, vis, conf, False)
assert_(_l2error(res0, res1) < epsilon)
......@@ -234,7 +251,7 @@ def test_against_dft(nxdirty, nydirty, epsilon, nchan, nrow):
@pmp("nrow", (10, 100, 1000))
@pmp("nchan", (1, 10))
@pmp("fov", (1,))
def test_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
def test_wstack_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
epsilon = 1e-7
np.random.seed(40)
pixsize = fov*np.pi/180/nxdirty
......@@ -252,8 +269,7 @@ def test_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
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))
res0 = np.zeros((nxdirty, nydirty), dtype=np.complex128)
mi, ma = np.min(uvw[:, 2]), np.max(uvw[:, 2])
nplanes = 10000
ws = mi + np.arange(nplanes)*(ma-mi)/(nplanes-1)
......@@ -264,29 +280,17 @@ def test_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
if len(iind) == 0:
continue
dd = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, iind, bl.ms2vis(ms, iind)))
res0 += conf.apply_wscreen(dd, 0.5*(ws[ii+1]+ws[ii]), adjoint=False).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)
res0 += conf.apply_wscreen(dd, 0.5*(ws[ii+1]+ws[ii]), adjoint=False)
res1 = _dft(uvw, vis, conf, True)
assert_(_l2error(res0, res1) < 1e-4) # Very inaccurate
n=np.sqrt(1-x**2-y**2)
nm1=n-1
for ii in idx:
phase = x*uvw[ii, 0] + y*uvw[ii, 1] + uvw[ii, 2]*nm1
res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
res1 /= n
assert_(_l2error(res0, res1) < 1e-4)
@pmp('nxdirty', [16, 64])
@pmp('nydirty', [64])
@pmp("nrow", (10, 100, 1000))
@pmp("nchan", (1,))
@pmp("fov", (50,))
def test_against_wdft_exact(nxdirty, nydirty, nchan, nrow, fov):
def test_wplane_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
epsilon = 1e-10
np.random.seed(40)
pixsize = fov*np.pi/180/nxdirty
......@@ -298,33 +302,44 @@ def test_against_wdft_exact(nxdirty, nydirty, nchan, nrow, fov):
speedoflight, f0 = 3e8, 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight)
# use exactly one random w
uvw[:,2] = uvw[0,2]
uvw[:, 2] = uvw[0, 2] # use exactly one random w
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)
mi, ma = np.min(uvw[:, 2]), np.max(uvw[:, 2])
nplanes = 1
w = mi
w = np.min(uvw[:, 2])
iind = ng.getIndices(bl, conf, flags)
dd = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, iind, bl.ms2vis(ms, iind)))
res0 = conf.apply_wscreen(dd, w, adjoint=False).real
res0 = conf.apply_wscreen(dd, w, adjoint=False)
res1 = _dft(uvw, vis, conf, True)
assert_(_l2error(res0, res1) < epsilon)
# 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.sqrt(1-x**2-y**2)
nm1 = n-1
for ii in idx:
phase = x*uvw[ii, 0] + y*uvw[ii, 1] + uvw[ii, 2]*nm1
res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
res1 /= n
assert_(_l2error(res0, res1) < 1e-4)
@pmp('nxdirty', [16, 64])
@pmp('nydirty', [64])
@pmp("nrow", (10, 100, 1000))
@pmp("nchan", (1, 10))
@pmp("fov", (1, 10))
@pmp("epsilon", (1e-2, 1e-7, 2e-12))
def test_wgridder_against_wdft(nxdirty, nydirty, nchan, nrow, fov, epsilon):
np.random.seed(40)
pixsize = fov*np.pi/180/nxdirty
conf = ng.GridderConfig(nxdirty=nxdirty,
nydirty=nydirty,
epsilon=epsilon,
pixsize_x=pixsize,
pixsize_y=0.5*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 = ng.vis2dirty_wstack(bl, conf, idx, vis)
res1 = _dft(uvw, vis, conf, True).real
assert_(_l2error(res0, res1) < epsilon)
Markdown is supported
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