Commit c27da4fe authored by Martin Reinecke's avatar Martin Reinecke
Browse files

add very simple demo

parent 8f44eb9b
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
# number of rows in the measurement set
nrow = 10000
# number of channels
nchan = 100
# Invent mock UVW data. For every row, there is one UVW triple
uvw = np.random.rand(nrow,3)
# scaling factor to obtain the UV of a row for a given channel from the UV
# in the measurement set.
# (This is proportional to 1/lambda, and future versions of the gridder will
# probably simply take the lambda array to make things more intuitive).
scaling = 1. + np.arange(nchan)*0.01 # just use silly values for this example
# Build Baselines object from the geometrical information
baselines = ng.Baselines(coord=uvw, scaling=scaling)
# Build GridderConfig object describing how gridding should be done
# urange and vrange are some very unintuitive quantities describing the extent
# of the UV grid. This will be replaced by something better; most likely we
# will ask for FOV of the dirty image and compute everything from that.
gconf = ng.GridderConfig(nxdirty=1024, nydirty=1024, epsilon=1e-7, urange=1., vrange=1.)
# 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, chbegin=3, chend=6, wmin=-0.1, wmax=0.1)
# 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()
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