Commit 6a41bee0 authored by Philipp Arras's avatar Philipp Arras
Browse files

Use scipy sparse matrix

parent fa6e1364
......@@ -345,64 +345,39 @@ def test_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
res1 /= n
assert_(_l2error(res0, res1) < 1e-4)
@pmp('nxdirty', (32,16))
@pmp('nydirty', (32,64))
@pmp("nrow", (100,1000))
@pmp("nchan", (1,5))
@pmp('epsilon', (1e-3,1e-6))
@pmp('nxdirty', (32, 16))
@pmp('nydirty', (32, 64))
@pmp("nrow", (100, 1000))
@pmp("nchan", (1, 5))
@pmp('epsilon', (1e-3, 1e-6))
def test_holo_from_correlations(nxdirty, nydirty, nchan, nrow, epsilon):
np.random.seed(420)
uvw = np.random.randn(nrow, 3)
uvw[:, 2] = 0
lightspeed, f0 = 3e8, 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
fmax = freq.max()
# Nyquist sampling
uv_max = np.abs(uvw * fmax/lightspeed).max()
cell = 0.75/(2*uv_max)
conf = ng.GridderConfig(nxdirty=nxdirty,
nydirty=nydirty,
epsilon=epsilon,
pixsize_x=cell,
pixsize_y=cell)
bl = ng.Baselines(coord=uvw, freq=freq)
flags = np.zeros((nrow, nchan), dtype=np.bool)
idx = ng.getIndices(bl, conf, flags)
# random vector to apply to
tmp = (np.random.randn(2*nxdirty, 2*nydirty) +
1.0j*np.random.randn(2*nxdirty, 2*nydirty))
# result of implicit operator
res1 = ng.apply_holo(bl, conf, idx, tmp)
# get explicit operator
from scipy.sparse import coo_matrix
np.random.seed(42)
bl, conf, idx = _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow)
nx, ny, gridsize = 2*nxdirty, 2*nydirty, 4*nxdirty*nydirty
grid = np.random.randn(nx, ny) + 1.0j*np.random.randn(nx, ny)
res1 = ng.apply_holo(bl, conf, idx, grid)
W = conf.W()
holo_d = {}
for i in range(0, W):
for j in range(0, W):
# add entry to dict
holo_d[(i,j)] = np.fft.fftshift(ng.get_correlations(bl, conf, idx, du=i, dv=j)).flatten()
if i:
holo_d[(i,-j)] = np.fft.fftshift(ng.get_correlations(bl, conf, idx, du=i, dv=-j)).flatten()
# apply to random vector
Ntot = 2*nxdirty * 2*nydirty
x = np.fft.fftshift(tmp).flatten()
res2 = np.zeros(Ntot, dtype=np.complex128)
for key in holo_d:
i, j = key
shift = i*2*nydirty + j
diag = holo_d[key]
if shift:
# upper diagonal
res2[0:Ntot - shift] += diag[0:Ntot - shift] * x[shift::]
# lower diagonal
res2[shift::] += diag[0:Ntot - shift] * x[0:Ntot - shift]
else:
res2 += diag * x
res2 = np.fft.ifftshift(res2.reshape(2*nxdirty, 2*nydirty))
assert_array_almost_equal(res1, res2, decimal=8)
\ No newline at end of file
mat = np.zeros(2*(gridsize,))
n = gridsize*(2*W-1)**2
row = np.zeros((n,))
col = np.zeros((n,))
data = np.zeros((n,))
ind = 0
foo = np.arange(0, gridsize)
for ii in range(-W+1, W):
for jj in range(-W+1, W):
ilow = ind*gridsize
ihigh = (ind+1)*gridsize
data[ilow:ihigh] = ng.get_correlations(bl, conf, idx, du=ii, dv=jj).ravel()
row[ilow:ihigh] = foo
ycoord = foo % ny
xcoord = foo-ycoord
tmp = (ycoord+jj) % ny + xcoord
ycoord = tmp % ny
xcoord = tmp-ycoord
col[ilow:ihigh] = (xcoord + ii*ny) % gridsize + ycoord
ind += 1
mat = coo_matrix((data, (row, col)), shape=(gridsize, gridsize))
res2 = mat.dot(grid.ravel()).reshape(nx, ny)
assert_allclose(res1, res2)
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