Commit 87d0487a authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'simpler_arrays' into 'master'

Merge new developments into master

See merge request !4
parents 7b579dbf e631d7e8
import numpy as np
import nifty_gridder as ng
# Some generic nomenclature:
# (If this sounds completely stupid, it's because I have no radio background
# whatsoever; I'm more than happy to change this where needed!)
#
# nrow : integer (number of rows in a measurement set)
# nchan : integer (number of channels in a measurement set)
# baselines: an object containing
# - a float(nrow,3) array; this are the UvW coordinates in the measurement set
# - a float(nchan) array; this describes how the UVW need to be scaled for every channel
# gconf: an object containing information about a gridding setup,
# i.e. resolution of dirty image, requested accuracy etc.
# ms: a complex(nrow, nchan) array containing the visibilities of a measurement set
# Later on this will probably become (nrow, nchan, npol)
# flags: a bool(nrow, nchan) array. Where True, the corresponding visibilities
# will be ignored.
# idx: a 1D integer array containing indices of selected visibilities. One index
# ranges from 0 to nrow*nchan-1 and encodes row and channel number simultaneously
# to save space.
# vis: a 1D complex array, which is always accompanied by an "idx" array. It contains
# the visibilities of a "ms", extracted at "idx".
# grid: oversampled 2D grid in UV space onto which the visibilities are gridded
# (resolution is higher than that of the dirty image)
# dirty: float(nxdirty, nydirty): the dirty image
f0 = 1e9 # rough observation frequency
npixdirty = 1024
pixsize = np.pi/180/60/npixdirty # assume 1 arcmin FOV
speedoflight = 3e8
# number of rows in the measurement set
nrow = 10000
# number of channels
nchan = 100
# Frequency for all channels in the measurement set [Hz].
freq = f0 + np.arange(nchan)*(f0/nchan) # just use silly values for this example
# Invent mock UVW data. For every row, there is one UVW triple
# NOTE: I don't know how to set the w values properly, so I use zeros.
#uvmax = pixsize*np.max(freq)/sp
uvw = (np.random.rand(nrow,3)-0.5) / (pixsize*f0/speedoflight)
uvw[:,2] = 0.
# Build Baselines object from the geometrical information
baselines = ng.Baselines(coord=uvw, freq=freq)
# Build GridderConfig object describing how gridding should be done
# pixsize_x and pixsize_y are given in radians
gconf = ng.GridderConfig(nxdirty=npixdirty, nydirty=npixdirty, epsilon=1e-7, pixsize_x=pixsize, pixsize_y=pixsize)
# At this point everything about the experimental setup is known and explained
# to the gridder. The only thing still missing is actual data, i.e. visibilities
# and flags.
# Invent mock flags. This is a bool array of shape (nrow, nchan).
# For this test we set it completely to False
flags = np.zeros((nrow, nchan), dtype = np.bool)
# extract the indices for the subset of channels and w that we want to grid.
# The gconf object is needed here because knowing the gridding parameters helps
# to optimize the ordering of the returned indices.
# For parallel processing it is possible to create multiple index sets, each
# covering a different range of channels, or to generate all w slices in
# parallel.
# If the complete "flags" array does not fit into memory, we can adjust the
# interface: for example, we could just pass the flag sub-array that matches the
# selected channel range.
idx = ng.getIndices(baselines, gconf, flags)
# Invent mock visibilities. Currently we have just NPOL=1, resulting in
# (nrow, nchan) complex visibility values
# Currently, the gridder code refers to this as "ms". Suggestions with more
# appropriate names are welcome!
ms = np.random.rand(nrow,nchan)-0.5 + 1j*(np.random.rand(nrow,nchan)-0.5)
# extract the visibility data at the obtained indices from ms.
# For large-scale datasets where ms does not fit into memory, this needs to be
# done differently, but should still be straightforward.
vis = baselines.ms2vis(ms, idx)
# perform the gridding
# Many of these operations can be called in parallel with identical baselines
# and gconf arguments; they won't interfere with each other.
grid = ng.vis2grid(baselines, gconf, idx, vis)
# convert gridded data in UV space to dirty image (i.e. FFT, cropping,
# multiplication with correction function)
dirty = gconf.grid2dirty(grid)
# Adjointness test
dirty2 = np.random.rand(*dirty.shape)
ms2 = baselines.vis2ms(ng.grid2vis(baselines, gconf, idx, gconf.dirty2grid(dirty2)), idx)
print (np.vdot(ms,ms2).real, np.vdot(dirty, dirty2))
import matplotlib.pyplot as plt
plt.imshow(dirty)
plt.show()
This diff is collapsed.
This diff is collapsed.
import nifty_gridder as ng
import numpy as np
import pytest
from numpy.testing import assert_, assert_allclose
pmp = pytest.mark.parametrize
@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(nxdirty, nydirty, nrow, nchan, epsilon):
f0 = 1e9 # rough observation frequency
pixsize = np.pi/180/60/nxdirty # assume 1 arcmin FOV
speedoflight = 3e8
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)
gconf = ng.GridderConfig(nxdirty=nxdirty, nydirty=nydirty,
epsilon=epsilon, pixsize_x=pixsize, pixsize_y=pixsize)
flags = np.zeros((nrow, nchan), dtype = np.bool)
idx = ng.getIndices(baselines, gconf, flags)
ms = np.random.rand(nrow,nchan)-0.5 + 1j*(np.random.rand(nrow,nchan)-0.5)
vis = baselines.ms2vis(ms, idx)
grid = ng.vis2grid(baselines, gconf, idx, vis)
dirty = gconf.grid2dirty(grid)
dirty2 = np.random.rand(*dirty.shape)-0.5
ms2 = baselines.vis2ms(ng.grid2vis(baselines, gconf, idx,
gconf.dirty2grid(dirty2)), idx)
assert_allclose(np.vdot(ms,ms2).real, np.vdot(dirty, dirty2))
def test_pickling():
try:
import cPickle as pickle
except ImportError:
import pickle
import nifty_gridder as ng
# Have to use cPickle and pickle protocol 2 for this to work
# unfortunately
pickle_protocol = 2
gc = ng.GridderConfig(1024, 1024, 2e-13, 2.0, 2.0)
pickled = pickle.dumps(gc, pickle_protocol)
gc2 = pickle.loads(pickled)
assert_(gc is not gc2)
assert_(gc.Nxdirty() == gc2.Nxdirty())
assert_(gc.Nydirty() == gc2.Nydirty())
assert_(gc.Epsilon() == gc2.Epsilon())
assert_(gc.Pixsize_x() == gc2.Pixsize_x())
assert_(gc.Pixsize_y() == gc2.Pixsize_y())
try:
import cPickle as pickle
except ImportError:
import pickle
import nifty_gridder as ng
# Have to use cPickle and pickle protocol 2 for this to work
# unfortunately
pickle_protocol = 2
gc = ng.GridderConfig(1024, 1024, 2e-13, 2.0, 2.0)
pickled = pickle.dumps(gc, pickle_protocol)
gc2 = pickle.loads(pickled)
assert gc is not gc2
assert gc.Nxdirty() == gc2.Nxdirty()
assert gc.Nydirty() == gc2.Nydirty()
assert gc.Epsilon() == gc2.Epsilon()
assert gc.Pixsize_x() == gc2.Pixsize_x()
assert gc.Pixsize_y() == gc2.Pixsize_y()
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