Commit 1f7224ef authored by Martin Reinecke's avatar Martin Reinecke
Browse files

update README; remove old files

parent 11a10e5f
......@@ -7,12 +7,10 @@ Library for high-accuracy gridding/degridding of radio interferometry datasets
Programming aspects
- written in C++11, fully portable
- shared-memory parallelization via OpenMP
- shared-memory parallelization via OpenMP and C++ threads.
- Python interface available
- kernel computation is performed on the fly, avoiding inaccuracies
due to table lookup and reducing overall memory bandwidth
- gridding/degridding of a single w-plane is OpenMP-parallelized with good
scaling behaviour
Numerical aspects
- uses the analytical gridding kernel presented in
......@@ -22,3 +20,11 @@ Numerical aspects
- in combination these two aspects allow extremely accurate gridding/degridding
operations (L2 error compared to explicit DFTs can go below 1e-12) with
reasonable resource consumption
Installation and prerequisites
- Clone the repository
- execute `pip3 install --user .`. This requires the g++ compiler.
`numpy` and `pybind11` should be installed automatically if necessary.
For the unit tests, `pytest` is required.
The file `demo_wstack_realdata.py` requires `casacore` and measurement set
data.
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()
import matplotlib.pyplot as plt
import nifty_gridder as ng
import numpy as np
np.random.seed(40)
ndirty, nrow = 512, 6 # breaks
# ndirty, nrow = 64, 1 # almost there
nchan, fov, epsilon = 1, 10, 1e-7
pixsize = fov*np.pi/180/ndirty
conf = ng.GridderConfig(nxdirty=ndirty,
nydirty=ndirty,
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)
res1 = ng.vis2dirty_wstack(bl, conf, idx, vis).real
x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in [ndirty, ndirty]],
indexing='ij')
x *= conf.Pixsize_x()
y *= conf.Pixsize_y()
res0 = np.zeros((ndirty, ndirty))
eps = np.sqrt(x**2 + y**2)
s = np.sin(eps)
nm1 = -s*s/(1 + np.cos(eps))
n = nm1 + 1
for ii in idx:
phase = x*uvw[ii, 0] + y*uvw[ii, 1] + uvw[ii, 2]*nm1
res0 += (vis[ii]*np.exp(2j*np.pi*phase)).real
res0 /= n
plt.imshow(res0)
plt.colorbar()
plt.savefig('groundtruth.png')
plt.close()
plt.imshow(res1)
plt.colorbar()
plt.savefig('new.png')
plt.close()
plt.imshow(res1/res0)
plt.colorbar()
plt.savefig('ratio.png')
plt.close()
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