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

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