3. Imager minor-cycle refactoring

3.1 'Clean' refactoring plan for casa-3.5

Parent JIRA ticket.

Approximate Class layout : 1

Backward-compatibility : Make new code available as a new task 'tclean'. Once it is completely ready, both tasks will exist for a while, while scripts are being converted over. Then tasks switch to 'oldclean', 'clean' and stay for a little longer before 'oldclean' goes away completely.

3.1.1 ImageSkyModelControl Class

Implement an ImageSkyModelControl class to handle iteration-control and calls to major/minor-cycle code. This is the part that is currently replicated in more than 6 XXXImageSkyModel classes, behaves inconsistently in places, and prevents easy plug/play of minor-cycle algorithms.

1. Main loops that control major-cycle calls : makeResiduals and makePSFs. This includes special cases such as niter=0 makes only residuals and no psfs, save model-to-ms in the first and last iterations and maybe in-between, don't recompute residuals if the user supplies a starting residual image, trigger Hogbom and Clark as special-cases of CS-clean where only the initial major cycle is done.

2. Iterations and stopping thresholds : niter (major and minor), user-specified threshold, minor-cycle-stopping-threshold (cyclefactor,psf-sidelobe,max-psf-fraction,peak-residual). This has to work correctly in interactive mode too, where the user can change niter and user-threshold.

3. Iterations on fields, stokes, channels : Separate iteration-counters, rules for stopping-thresholds (per plane or combined), controls to prevent unnecessary model-predictions when there are no components to predict.

4. Mask behaviour : Rules for merging across multiple fields, extending masks across stokes, channels, plugin for autoboxing.

5. PSF validation : Currently, problems at most data-selection stages results in 'bad PSF' exceptions from this part of the code, when trying to fit a clean-beam. Add an explicit PSF validation with helpful error-messages, and a clean-return.

6. Runtime parameter modification : Ability to modify iteration parameters while clean is running, using the AppRC class that operates via an editable text file. This is in addition to what the interactive-clean GUI provides at the end of every 'npercycle' iterations. (Eventually, it would be good to have the GUI stay alive and enable runtime param modification through it).

7. Log messages : Iteration counters and thresholds per field/stokes/channel, interactive and non-interactive, flux gathered in current set of minor-cycles and total, etc...

CSSC Tests : Test the above (only iteration-control and messages) via one of the following : (a) a new option for one of the existing task parameters (perhaps 'imagermode'/'psfmode'). If possible, we can temporarily connect one of the minor-cycle cleaners (ms-clean), so that it can be tested with images. (b) first implement a light-weight task interface following the new choices, and test it directly there, as a separate 'tclean' task.

3.1.2 Deconvolver class

Merge minor-cycle cleaners (used by task_clean/Imager) with the existing deconvolver class (used by task_pclean/Deconvolver) so that they all use the same code. Currently, only ms-clean is shared by both Imager and Deconvolver.

1. Merge minor-cycle cleaners : Add cleaners for hogbom, clark and ms-mfs to the Deconvolver. This will make it usable from task_pclean. Make Imager use Deconvolver for all minor-cycle algorithms.

2. Synchronize iteration-control and threshold-checks for all deconvolvers : This is mainly to give synchronized log messages for 'amount of flux found', etc. Preferably in a Deconvolver base class.

3. Image-restoration options : Add a function to the Deconvolver base class for default 'restoration', to replace the code in Imager2.cc. Any derived minor-cycle algorithms can do their own types of 'restoration' if needed (e.g. ms-mfs makes spectral-index maps too). Cube cleaning can get channel-dependent restoring beams.

4. Mask modification options : Derived deconvolvers can treat input masks differently. For example, multi-scale cleaners expand masks per scale. Wide-band masks could be merged across terms, etc...

5. Memory-counter : Factor out all memory-allocation calls into one function per deconvolver. This is for us to keep track of how much mem-usage we expect. Also include a memory-counter function in ImageSkyModel - the class that holds all the lists of images that get passed around everywhere.

6. Runtime parameter-modification : Place checks on the AppRC state within the cleaners at the appropriate place. Preferably in a Deconvolver base class.

CSSC Tests : Test all the minor-cycle algorithms. This should be possible via the current clean-task interface. If not, some new parameters will be added temporarily.

3.1.3 New clean task-interface

Make it more logical than the current interface. Reconcile imagermode, psfmode, gridmode, etc. : See section . We may not get to this for the 3.5 release.

1. Empty task : An empty clean task, with all the parameters, and messages for what all combinations will do. This should handle multi-field, multi-term, diff ftmachines, single-plane, cubes, ability to run parts in isolation (make dirty image, make pb, etc..). This will also sync with the requirements of task_pclean, so that we can eventually have only one main task.

2. Connect new task to ImageSkyModelControl : Use the existing tool-interface and Imager.cc to connect the new task to ImageSkyModelControl. Re-use the existing Imager.cc and make in-place changes wherever possible, otherwise make copies of functions and delete the old ones when we move over. Try to re-use python task code (i.e. cleanhelper), but wherever possible, re-implement it in C++.

3. Log messages : Imager-setup messages about data-selection (number of visibilities, frequency-range,etc), image-parameters (freq/velo ranges, etc), expected image rms, etc....

4. Connect to FTMachines : Most of this already exists/works (CubeSkyEquation class). Enable all sensible combinations of minor-cycle deconvolvers with FTMachines (ft, mosaic, wproj). Primary-beam correction will be restricted to post-deconvolution image-domain correction in python at the end of 'clean'.

5. Write 'unit'-tests : All ftmachine/devonvolver combinations, iterations for fields,stokes,channels, etc....

CSSC Tests : Task interface, setup log-messages, test with all ftmachines. Also, some users have asked for small tasks that do only one thing , for example "make dirty image" or "make pb", etc. It would be best if CSSC/user-support writes small wrapper tasks, using the buildmytasks framework, and call 'clean' and 'cleanhelper' functions internally.

3.2 New Clean task parameters

This section lists an alternate tasking interface for image-reconstruction in CASA. This interface follows the underlying structure of the algorithms being used, and splits the parameters into five sections : data selection, image-definition, iteration-control, imaging, deconvolution.

It will be possible to run only 'selection + image-definition', generate some logger output ( Number of visibility points for image-rms calculations, selected frequency-ranges in the data-frame and LSRK, etc.. ) and then exit. This is to allow verification of image-coordinates before proceeding with the full imaging.

A run with niter=0 would do one major cycle to make PSFs and residual images (and return a dictionary of info). If a modelimage is specified, this step will save or predict the model (the equivalent of task 'ft' but with full flexibility of setup and ftmachine choices). Irrespective of 'niter', this first major cycle will generate logger information about sum-of-weights, expected image-rms, min and max spatial-scales spanned by the uv-coverage, clean-beam size, PSF sidelobe level, etc. If possible, if no parameters are changed, it should be possible to re-start deconvolution directly from the minor cycles, and avoid re-doing the major cycle.

3.2.1 Data Selection

Select a subset of a dataset (or multiple datasets). All selected data will be used for image reconstruction.

-> Imager tool interface : im.open(); im.selectvis();

vis

'' or ['','']

Name of input visibility file (for multiple ms, give a list of strings)

field

'' or ['','']

Field Name or id (for multiple ms, give a list of strings)

spw

'' or ['','']

Spectral windows e.g. '0$\sim$ 3', '' is all (for multiple ms, give a list of strings)

selectdata

True

More selection (for multiple ms, give a list of strings)

timerange	'' or ['','']	Range of time to select from data
uvrange	'' or ['','']	Select data within uvrange
antenna	'' or ['','']	Select data based on antenna/baseline
scan	'' or ['','']	Scan number range

False

No more data selection

Keep the 'selectdata' parameter even though its default is True and it doesn't help save screen-real-estate. Jeff has a plan to consolidate 'selection' across tasks, so until then, don't change anything.

3.2.2 Image Definition

Define the coordinate systems of output images.

-> Imager tool interface : im.defineimage();

Postpone discussion about mode='cube' vs mode='channel','velocity','frequency', as the CSSC does not yet agree on the 'cube' idea/document circulated by JO. For now, no changes here.

The following is a list of points to discuss, based on suggestions from KG and DP.

1. Remove 'outframe' because Imager always grids in LSRK, and provide an external tool/task to 'set' the conversion-layer in the image-coordinate system for use by the viewer and fits-writers for on-the-fly conversion.
2. Alternatively, add support for 'native-gridding' vs 'LSRK' gridding, to support datasets that are pre-gridded to the target image-channelization, and which do not need/want the LSRK conversion.
3. 'reffreq' is decided purely by data-selection (middle of the range in LSRK)
4. restfreq' : Needed before conversion to FITS. What about multiple rest-freqs ?
5. Perhaps remove support for velocity-image-coords, and 'veltype', and force this step as a post-deconvolution conversion ?
6. Add a new parameter 'excludefreqs' either here, or in data-selection. It would be a list of frequency ranges (or velocities), with units and ref-frames, to exclude from gridding/imaging. This is to prevent the user having to guess/calculate which TOPO data channels would map to these LSRK frequencies or velocity-ranges at at-least one timestep.

Postpone decisions about the existence of parameters 'reffreq', 'restfreq', 'outframe', until after KG's suggestions for reducing user-confusion about spw-selections vs image-freq-definitions, especially when frequency-frame changes require frequency-regridding, and there is NO one-to-one mapping between data channels and image channels.

mode

'mfs'

Combine all frequencies during imaging to produce a wideband image. The bandwidth is computed from the selected data-range, converted to LSRK.

reffreq	' '	Reference frequency of the continuum image. By default, it is the middle of the frequency-range computed from the data-range converted to LSRK. This is also the ref-freq used by the 'msmfs' deconvolution algorithm for its Taylor-polynomial model of the source spectrum.
nterms	2	Number of Taylor-polynomial terms used to model the sky spectrum.

'channel'

Cube, image-channels are defined via data-channel-frequencies converted to LSRK at first timestep.

nchan	-1	nchan
start	''	start chan (from zero in the MS, ignoring data-selection)
width	1	step size
interpolation	'linear'	Interpolation to use when regridding visibilities from data to image-frequencies. As a default, this should be 'nearest' for mode='channel', but we currently cannot have different subparam defaults for different parent-params
outframe	'LSRK'	The image is always made in LSRK. An external frame-convertor can be set to one of the following 'LSRK' 'LSRD' 'BARY' 'GEO' 'TOPO' 'GALACTO' 'LGROUP'.

'frequency'

Cube, with image-channels defined via frequency-ranges.

nchan	-1	nchan
start	''	start freq
width	1	step size
interpolation	'linear'	Interpolation to use when regridding visibilities from data to image-frequencies.
outframe	'LSRK'	The image is always made in LSRK. An external frame-convertor can be set to one of the following 'LSRK' 'LSRD' 'BARY' 'GEO' 'TOPO' 'GALACTO' 'LGROUP'.

'velocity'

Cube, with image-channels defined via velocity-ranges.

nchan	-1	nchan
start	''	start velocity
width	1	step size
interpolation	'linear'	Interpolation to use when regridding visibilities from data to image-frequencies.
outframe	'LSRK'	The image is always made in LSRK. An external frame-convertor can be set to one of the following 'LSRK' 'LSRD' 'BARY' 'GEO' 'TOPO' 'GALACTO' 'LGROUP'.
veltype	'radio'	Velocity definition of output image. 'optical' or 'radio'.

stokes

'I'

Stokes params or correlation products to image : I,Q,U,V, RR, LL, XX,YY and all combinations allowed by the input data (I,Q,U,V,IQ,IV,QU,UV,IQU,IUV,IQUV,RR,LL,RRLL,XX,YY,XXYY).

imagename

' '

Pre-name of output images ( [ ' ', ' ' ] for multi-fields )

imsize

[256, 256]

x and y image size in pixels. Single value: same for both. ( [ [,], [,] ] for multi-fields )

cellsize

['1.0arcsec', '1.0arcsec']

x and y cell size(s). Default unit arcsec. ( [ [' ',' '],[' ',','] ] for multi-fields)

phasecenter

' '

Image center: direction or field index ( [ ' ', ' '] for multi-fields)

mask

' ' or [ , , ]

Cleanbox(es), mask image(s), region(s), or a PB-level ( [ [ , , ] , [ , , ] ] for multi-fields)

modelimage

' ' or [' ' , ' ']

Name of model image(s) to initialize cleaning ( [ [ , , ] , [ , , ] ] for multi-fields)

outlierfile

' '

Outlier file (new format) with imagename, imsize, cell, phasecenter, mask, modelimage. A single main field is specified in the task-parameters, and only outliers are in this file. Outlier files cannot be mixed with list-modes.

(1) A convertor can be provided for the AIPS format, but without 'modelimage' support. i.e. no mixing of inline-lists and outlier-files for parameters un-supported by the AIPS format. If the convertor is run separately, the user can then edit the files to add model-image info.

(2) Old CASA-Imager outlier-file format will go away.

Example outlierfile format : Two outliers, supplying a mask and modelimage only for the first. These two fields are in 'addition' to the main single field specified in the task parameters.

    imagename = 'M1_1' 
    imsize = [128,128]
    phasecenter = 'J2000 13h30m52.159 43d23m08.02'
    mask = ['out1.mask', 'circle[[40pix,40pix],5pix]' ]
    modelimage = 'out1.model'

    imagename = 'M1_2' 
    imsize = [128,128]
    phasecenter = 'J2000 13h24m08.16 43d09m48.0'

Inline lists for multi-field specification : Allow only for simple cases, but force outlier files as the only option for modes that require multiple levels of nested lists, which can also be different levels for different parameters (for example, multi-term + multi-field + multiple-masks per field).

Initial implementation will allow only outlier files and no inline-lists. Support for inline lists will then be added for 'specific, simple' cases for ease-of-use.

3.2.3 Iteration Control

Choose the type of iterations (major and minor cycles), gain, and thresholds. All iteration-counts are maintained per-field, per-channel and per-stokes.

-> im.clean()

niter

160

Maximum number of minor-cycle iterations, counted across major-cycles. This is over-ridden by the user-specified stopping-threshold.

npercycle

-1

Force hogbom-clean, with major-cycles only at the beginning and end. Do not use FluxLimit.

50

Maximum number of minor-cycle iterations before triggering a major-cycle. This is overridden by 'FluxLimit' calculations controlled by sub-parameters, to trigger a major cycle. Also overridden by the user-specified stopping-threshold.

Total number of major cycles $>=$ niter/npercycle

For example, niter=160, npercycle=50. Due to the flux-limit, an extra major-cycle is triggered after iteration 28. Major-cycles will occur after iterations 0, 28, 78, 128, 160.

FluxLimit = peakResidual x fractionOfPsf where fractionOfPsf = MIN(maxPeakFraction, MAX( minPeakFraction, PSFsidelobe x cyclefactor) )

This FluxLimit is controlled by the PSF sidelobe level and 'cyclefactor'

maxpeakfraction	0.8	The highest fraction of the peak residual that FluxLimit can be. Must be smaller than 1.0 !
minpeakfraction	0.1	The lowest fraction of the peak residual that FluxLimit can be. Must be larger than 0.0
cyclefactor	1.0	Scale-factor for the psf-sidelobe level. Rename to 'sidelobescalefactor' ?
nitertol	0.05	If we get to npercycle iterations, or a fluxlimit, and less than 5% of total iterations remain, stop and tell the user those last-few iterations won't be done. If set to 'zero', then those last few iterations will be done after the major cycle. May be dangerous to just do those extra iterations as it may go too deep.
checkdivergence	0.5	If the peak residual increases by this fraction, do not use the last model, but trigger a major cycle
negcomponent	-1	If larger than zero, stop minor cycle iterations when this many negative components are found. When 'multiscale' is turned on, this count applies to negative components at the largest spatial scale.

gain

0.1

Loop gain for cleaning

-1

Adaptive loop gain, following the Levenger-Macquart method, to go faster when possible, and slow down when in a shallow part of chi-square space (not discussed yet, but already available/tested for msmfs in current task)

threshold

'0.1mJy'

Main stopping threshold for each field, channel and stokes-plane. If any field, channel, or stokes-plane reaches this threshold, stop minor cycles for the current field, channel or stokes-plane, and move on to the next field, channel or stokes-plane. This overrides all other iteration-control.

interactive

True

Pop up a GUI for mask-drawing with (stop, continue-till-end, continue until next stopping-point). Stop after every major cycle. Also allow modification of niter, npercycle, maxpeakfraction, minpeakfraction, cyclefactor, threshold, gain. This will be implemented via a text-file on disk, but will have a GUI (for ease of use).

If ftmachine='refpb' or 'aproject', always overlay a contour of the PB at 'pblimit'. If pbcor = True, overlay a contour of the PB at 'pbmin', in a different colour.

If calready=True, save the model at every stopping-point, so that one can exit at any time.

pauseoption	False	Stop only at major-cycle boundaries
	True	A button on the interactive-clean GUI stays alive, and can be clicked to 'pause', even between minor-cycle iterations. Then, masks can be modified and all those parameters can also be modified via the GUI. An exit from with minor-cycle iterations, will trigger one last major cycle.
plotprogress	False	Nothing
	True	Show continuously-updating plots of total flux and peak residual

False

Run without a GUI, but allow runtime parameter modification via a text file on disk containing parameters. This will be the same text-file involved in parameter-modification via the GUI. File name will be [imagename].control. Decide later if 'plotprogress' should be possible here too.

To start with, do not implement `` chaniter : Clean each channel (minor and major cycles) completely before proceeding to the next ``. If later required (purely as a usage-mode and not for performance), it can be implemented at the top level, with the expected caveat of diminished performance, because of (a) ignoring data-locality in tiles, and (b) having to read multiple channels even to make one image channel because of freq-frame conversions.

An option to do a flat-sky normalization before showing the current residual image was requested. It is possible to do this extra normalization, but to avoid making an extra copy of the whole image, it would mean an extra sequence of divisions/multiplication. This is an extra step, and with small numbers near the edge, may lead to numerical instability. Therefore, the initial implementation will only 'show' PB-contours at the 'pblimit' level when ftmachine='refpb' and 'aproject'. The extra division/multiplication by the PB can be added later if still required.

3.2.4 Imaging Options

Choose weighting schemes and direction-dependent effects to apply while gridding visibilities and constructing the image to be sent to the deconvolution algorithm.

-> Imager tool interface : im.weight();im.setvp();im.setoptions();

Note : Mosaicing (phase-gradients on gridding-convolution functions) is implicit in all of the following, and is automatically triggered when multiple pointings are present in the selected data (i.e. no parameters required). To do a linear-mosaic with separate deconvolutions per facet, use multi-field options followed by im.linearmosaic()).

weighting

'natural'

Natural weighting

'uniform'

Uniform weighting

'briggs'

Robust weighting

robust	0.0	go between uniform and natural
npixels	0	npixels to decide uv-cell size over which weights are computed

'briggsabs'

Robust weighting - 2

robust	0.0	go between uniform and natural
noise	'1.0Jy'	noise estimate
npixels	0	npixels to decide uv-cell size over which weights are computed

uvtaper

False

No uv-tapering

True

Apply additional uv tapering of visibilities

outertaper	[]	uv-taper on outer baselines in uv-plane
innertaper	[]	uv-taper in center of uv-plane ?? (not implemented)

ftmachine

'ft'

Use the standard prolate-spheroidal for gridding.

'sd'

Use the FT of the single-dish PB as the gridding-convolution kernel

'refpb'

Use the FT of antenna-dependent PBs, at ref-freq, and one time as the gridding-convolution kernel current ftmachine=mosaic. Data weights are not considered in PB computation. Normalization for the minor-cycle is always flat-noise.

pblimit

0.02

limit beyond which to ignore

'aproject'

Use time-varying, frequency-dependent, and antenna-dependent gridding convolution functions constructed as FTs of PBs, rotated and scaled for time and freq variations, squinted, and weighted according to data-weights. Normalization for the minor cycle is always flat-noise..

pblimit	0.02	limit beyond which to ignore
cfcache	'sss'	cache directory

wprojplanes

1

No w-projection. Use the gridding convolution function selected via 'ftmachine'.

64

Use w-projection kernel (FT of Fresnel propagator) with this-many w-planes, in addition to whatever selected by 'ftmachine'

facets

1

For each image-field, use a set of 'facets x facets' sub-images during gridding. Deconvolution will always see a joint image (covering the area : 'facets x facets'). Faceting can be used as an alternative to w-projection (wprojplanes=1), along-with W-projection (wprojplanes$>$ 1), or along-with A-projection or 'ft'. In all cases (especially mosaicing), the convolution-functions used will calculated for the facet-size field-of-view.

pbcor

False

Output image represents $I^{sky} \times PB$ . No extra normalization.

True

Output image represents $I^{sky}$ . Do a post deconvolution correction using an average PB. If nterms$>$ 1, this does a post-deconvolution wide-band pb-correction.

pbmin

0.1

Mask for post-deconvolution pb-correction

Initially, not all possible valid combinations will be enabled. Un-supported combinations will be checked-for and stopped at the start of the task, and enabled when the internals are ready for it. The main point here is that the interface will support all physically-valid combinations in a way that follows the mathematics of the imaging process.

3.2.5 Deconvolution Options

Choose the deconvolution algorithm based on sky parameterization. This will run separately per field, channel, and stokes-planes (except when combinestokes=True). Iterations are the number of minor-cycle iterations decided by iteration-control. This part will also be available at the tool-level as stand-alone image-based deconvolvers.

-> Imager tool interface : im.setscales(), im.settaylorterms(), im.clean();

deconvolver

'hogbom'

A variant of the Hogbom clean algorithm. Use a point-source flux model. Use full PSF.

'clark'

A variant of the Clark-clean algorithm. Use a point-source flux model. Use truncated PSF, minor cycle includes periodic subtraction using gridded visibilities.

'clarkstokes'

Similar to 'clark', but search for peaks in $I^2 + Q^2 + U^2 + V^2$ . Requires a 4-stokes-plan input image.

'multiscale'

Model the sky as a collection of Gaussian-like components of different sizes.

scalevector	[0,3,10]	List of scale sizes in units of pixels, or as strings : ['2.0arcsec', '6.0arcsec','20.0arcsec'].
smallscalebias	0.6	Give more weight to small scales

'mem'

maximum entropy (SUBPARAMS)

'asp'

asp-clean (SUBPARAMS)

autobox

False

No autoboxing

True

Do autoboxing. Similar to interactive=True one could contol whether the autoboxing should be repeated every N minor-cycle iterations or at major-cycle boundaries. (Add autoboxing SUBPARAMS)

restoringbeam

'' or ['','','']

Gaussian restoring beam for CLEAN image. Default is to use a frequency-dependent beam, fitted from the PSF main lobe at each frequency.

calready

False

Do not save a model anywhere inside the MS.

True

Save a model for later use in self-calibration or plotting.

usescratch	False	Save model as a virtual model-column record. The ftmachine state, with the image and data-selection info are stored.
	True	Write to the model-data scratch column in the MS.

async

False

If true the taskname must be started using clean(...)