imcontsub¶

The purpose of this software is to subtract continuum in a parallel/distributed environment or on a single computer system. The software leverages MPI, however can be run on a simple laptop or a large supercomputer.

The imcontsub task was originally based on the version of the ACES robust_contsub.py script from Dec 2019. It is still available by specifying iterativeclip=false. It performs the following processing for each spectrum:

exclude outliers (top and bottom 5%),

do approximate linear fit of the spectrum (using every 10th data point)

subtract off the linear fit to remove source spectral index to 1st order

determine mask based on threshold and IQR relative to median (see below)

use mask to fit original channels with polynomial of given order

subtract fit and write out residual spectrum

In Oct 2023 the algorithm was improved to implement iterative robust clipping. The new algorithm consists of the following steps:

exclude outliers (top and bottom 5%),

do linear fit to remaining points

exclude the 50% largest residuals and redo linear fit

determine a mask based on threshold and IQR relative to median

(on the first pass use minimum of Q50-Q25 and Q75-Q50 to estimate IQR)

repeat above steps until there are no further changes (up to 5 times)

use mask to fit original channels with polynomial of given order

subtract fit and write out residual spectrum

Running the program¶

It can be run with the following command, where “config.in” is a file containing the configuration parameters described in the next section.

$ <MPI wrapper> imcontsub -c config.in

Note that the number of ranks needs to divide evenly into the 2nd axis of the input cube (Nx,Ny,Nz). Each rank processes a block of Ny/Nranks x-z planes (containing Nx * Ny/Nranks spectra of length Nz each).

Configuration Parameters¶

Parset parameters understood by imcontsub are given in the following table (all parameters must have imcontsub prefix, i.e. imcontsub.order).

Parameter	Type	Default	Description
inputfitscube	string	None	Image cube to work with. The input file has to be in FITS format. It is read and left unchanged.
outputfitscube	string	generated	The name of the cube to write the results to. If the file already exists, it will be overwritten. If the output name is left unspecified, it is generated from the input name: mycube.fits -> mycube.contsub.fits
order	int	2	The order of the polynomial used to fit for the continuum.
threshold	float	2.0	The threshold (in robust rms units) used to decide which channels to include in the fit for the continuum. A channel is included if it is not flagged and: abs(pixval-median) < threshold*rms where rms is estimated from the IQR (inter quartile range)
blocksize	int	0	Do the subtraction in blocks of channels, this serves two needs: match the beamforming interval to deal with discontinuities at boundaries; better model the continuum by fitting over a smaller frequency range. The default of 0 will fit and subtract the whole spectrum at once
shift	int	0	Shift the origin of the subtraction blocks ‘left’ by this many channels to match the location of beamforming steps.
interleave	bool	false	Interleave the fit/subtract blocks, using blocksize channels to fit the spectrum, but subtracting only the central 50%, stepping by half a block. This avoids some edge effects of the fitting process obvious in the rms spectrum
iterativeclip	bool	true	If true, use iterative robust clipping to reject emission and absorption features before fitting continuum. If false, use the previous algorithm which was less accurate and mainly aimed at rejecting emission features
imageacces	string	collective	The default is to use MPI collective IO, which is the most efficient way to read and write the data. If that is not available, specify ‘individual’ for all ranks to work independently and write in sequential mode.
imageHistory	vector<string>	none	Specify lines to add to the image history

Example¶

# The input fits cube
imcontsub.inputfitscube                           = mycube.fits
# The output fits cube
imcontsub.outputfitscube                          = mycube.contsub.fits
# The polynomial order - doing a linear fit
imcontsub.order                                   = 1
# The threshold - include channels within 2.5 sigma
imcontsub.threshold                               = 2.5
# The blocksize - match default beamforming interval
imcontsub.blocksize                               = 54