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 is based on the version of the ACES robust_contsub.py script
from Dec 2019.
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

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


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                                                  |
+------------------------+------------+------------+----------------------------------------------------------+


Example
-------

.. code-block:: bash

    # 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