We have tested our CUDA code on a variety of simulated datasets. Some examples from these simulations are shown below.

System requirements:

• CUDA enabled NVIDIA GPU
• pycuda
• scikits.cuda
• pyfits

Download the code:

• gICLEAN.py
• PDF readme document

We have tested our CUDA code on a variety of simulated datasets. Some examples from these simulations are available here.

Gridding:

Synthesis aperture arrays sample the Fourier transform of the sky brightness towards the source. In order to make a science image you must first invert these complex valued measurements, V(u,v), by taking a Fourier transform. The common technique is to grid the sampled data points and evaluate the discrete Fourier transform via an FFT. We grid using a convolution with a standard spheroidal function and use the CUDA SciKit wrapper for cuFFT to calculate the sky image and point source response function. We also include a rudimentary weighting scheme that can be used to accentuate different aspects of the data. We use the following kernels:

• gridVis_wBM_kernel: Each pixel in the uv plane goes through the data and check to see whether the pixel is included in the convolution. This kernel also calculates the point spread function and the local sampling from the data (for applying the weights later).
• wgtGrid_kernel: This will apply the weights.
• dblGrid_kernel: To improve speed, the previous kernel is only evaluated for half the grid. Since the data is Hermitian (the sky brightness is purely real) we can figure out what the other half looks like afterwards.
• nrmGrid_kernel/nrmBeam_kernel: This normalizes an array.
• corrGrid_kernel: The corrects for the convolution we applied when gridding.
• shiftGrid_kernel: This shifts a grid for computing an FFT.
• trimIm_kernel: We pad the uv grid to improve image quality and this kernel trims it to the requested image size.

Hogbom CLEAN algorithm:

Radio antenna arrays typically have non-uniform spacing, which gives them access to low and high frequencies simultaneously, but with sparse sampling in the UV plane. A Fourier tranform of this sparsely sampled data yields significant artifacts in the resulting image. For example, an image of a point source will have many sidelobes throughout the 2D image, which can reach intensities within an order of magnitude of the primary lobe.

This point source response pattern is called the dirty beam. The mage produced by convolving an astronomical source with the dirty beam is called the dirty image or dirty map. A clean beam is typically constructed by fitting a 2D Gaussian to the primary lobe of the dirty beam. If the locations of all astrophysical sources in the dirty image are known, a clean image can be approximated by convolving this model with the clean beam.

Hogbom (1974, Astronomy and Astrophysics, 15, 417) presented a classic image deconvolution algorithm used in radio astronomy to remove these instrumental artifacts. We implement the Hogbom algorithm on the GPU using pyCUDA (“cuda_hogbom” function). The algorithm uses 3 CUDA kernels:

• find_max: The Hogbom algoritm begins by finding the highest intensity point (in terms of absolute value) on the dirty image. Our implementation uses the gpuarray.max an atomic exchange operator to find the maximum, and then the find_max kernel is called to identify the array index of the maximum and record the maximum value at this position on the image model.
• sub_beam: In each iteration, the dirty beam is scaled, positioned at the maximum point found in step 1, and subtracted from the dirty image. Our sub_beam kernel performs this subtraction on the gpu, and simultaneously adds the scaled and shifted clean beam to the clean image.
• add_noise: The final step in the Hogbom algorithm is to add the residuals from the dirty image, after iterative beam subtraction, to the model cleaned image as ‘noise.’ We do this using an element wise pyCUDA kernel. Data input

gICLEAN.py is executed as a python script from the command line:

python gICLEAN.py [example] [image size] [make figures]

The command line options can also be edited within the main routine:

• [example]: Choose a dataset to image. We provide the options of a ‘gaussian’, ‘ring’, ‘mouse’, or ‘hd163296’. All except the last are simulated datasets. The default is ‘gaussian’.
• [image size]: The number of pixels on an image side. The default is 1024.
• [make figures]: 0/1 toggles whether or not summary images are created. The default is 0=False.

This will compile and run the CUDA kernels and generate final images. It requires an ALMA style visibility dataset in FITS format that is set by the vfile keyword. This version is good for single channel images only (although it is in principle scalable for spectral lines or multi-frequency synthesis). We have supplied several examples along with reasonable imaging parameters in our code with the data available in the examples section. The gaussian, ring, and mickey examples are noiseless, simulated datasets while hd163296 is a single channel of the CO J=3-2 line observed by the sub-millimeter telescope ALMA.

First, set up the environment (this applicable to users of the Harvard SEAS Resonance cluster only):

gpu-login

To perform gridding and CLEANing for the Gaussian source example, simply run:

python gICLEAN.py

This script takes a few command line options:

python gICLEAN.py [EXAMPLE] [ISIZE] [PLOTME]

• example: Use this to select which example dataset to operate on. Options are “gaussian,” “ring,” “mouse,” and “hd163296.” FITS files containing visibilities (raw data) for each example to be gridded and CLEANed are available in a tar archive in the examples section .
• ISIZE: The output image size, in pixels. This is ‘1024’ by default; to make e.g. an image four times larger on each size, choose ‘4096’.
• PLOTME: Turn this on to generate various png plots using matplotlib. This slows down the code, because it needs to grab data from the GPU and do the plotting, but is useful for debugging and data display.

Within this module, there are several functions that have their own configurable options, although these are not given command line switches.

cuda_gridvis parameters:

• settings: This is a Python dictionary containing the following imaging parameters:
• vfile: A string with the location of the dataset (see Data Input section)
• cell: The size of a single pixel in arcseconds.
• imsize: The number of pixels on an edge. Images are square and must be a power of 2.
• briggs: This controls the visibility weights. A larger number corresponds to better sensitivity while a lower number will improve resolution.
• plan: This is the plan for cuFFT to follow.

cuda_hogbom parameters:

• thresh: This defines the stopping threshold for the Hogbom algoirthm, Iterations continue until the maximum intensity point on the dirty image is a factor less than thresh of the maximum in the original dirty image. The default value is thresh=0.2.
• gain: This defines the “strength” of each iteration. At each iteration, the dirty beam is multiplied by gain before subtraction from the ditty image. By default, gain=0.1, meaning that 10 iterations are required to fully clean an image with a single point source.