UVMultiFit’s Documentation

A CASA-based flexible visibility-fitting engine developed at the Nordic node of the ALMA Regional Center.

Current development is partially founded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 730562 [RadioNet], in the frame of the RINGS Working Package (Radio Interferometry New Generation Software).

This program can be used to fit multi-component models to the visibilities in one (or several) measurement set(s). The fits can be performed to the continuum or in spectral-line mode. Advanced model structures can be implemented easily and several different fits can be performed with no need of re-loading the data.

Installation

Just follow these steps:

  1. Download the latest version from Launchpad.

  1. Untar the contents of the package into a directory of your choice (e.g., in your .casa directory). For example: 

    mkdir ~/.casa/NORDIC_TOOLS
mv UVMULTIFIT_*.tar.gz   ~/.casa/NORDIC_TOOLS
cd ~/.casa/NORDIC_TOOLS
mkdir UVMULTIFIT.current
tar -xf UVMULTIFIT_*.tar.gz -C ./UVMULTIFIT.current
The names /NORDIC_TOOLS/UVMULTIFIT.current are just a suggestion.

  1. Install the dependencies (i.e., GNU Scientific Library, so far).

  • In GNU/Linux, execute one of these commands: 

    sudo apt-get install libgsl2

sudo apt-get install libgsl-dev

sudo apt-get install libgsl0ldbl

  • In MacOS, execute all these commands: 

    sudo port install gsl
export LIBRARY_PATH="/opt/local/lib"
export LD_LIBRARY_PATH="/opt/local/lib"

  1. From the same terminal where you ran the previous step, go to the installation directory (if you are not already there) and compile the program: 

    cd ~/.casa/NORDIC_TOOLS/UVMULTIFIT.current
$CASABASE/bin/python setup.py build_ext --inplace

    After this step, the file _uvmultimodel.so should have been created.

    Note

    Check that this file has indeed been created by you, and is not the pre-compiled file coming with the package (if there is one). Executing ls -ltr should list it at the end and give you its exact creation date.

     


  1. Import the module into CASA. You can automate it by editing the file called init.py, which resides in your .casa directory: 

    gedit ~/.casa/init.py

    Add these two lines in your init.py: 

    UVMULTIFIT_PATH = YOUR_HOME_DIRECTORY/.casa/NORDIC_TOOLS/UVMULTIFIT.current
import imp
uvm = imp.load_source('uvmultifit', UVMULTIFIT_PATH+'/uvmultifit.py')

    where YOUR_HOME_DIRECTORY should be (guess it!) your home directory (i.e., the directory that is printed when you execute cd ~ ; pwd in a terminal).                    

 


  1. Note

    If there is a CASA error when loading the task, your local numpy version may be incompatible with the one from CASA. You may need to retrieve the numpy CASA version (given by np.__version__) and install it in your computer with the pip command (pip install numpy==CASA-VERSION-OF-NUMPY).

     


  1. Enjoy!

Any feedback, problem with the installation/running and/or bug report should be sent either to the ARC Nordic Node (contact@nordic-alma.se) or to the source maintainer (ivan.marti-vidal@chalmers.se).

Basic Usage

THESE ARE ALL THE KEYWORDS (AND THEIR DEFAULTS) OF THE CURRENT UVMULTIFIT VERSION: 

>>> uvmultifit(self,vis='', spw='0', column = 'data', field = 0, scans = [],
>>>          uniform=False, chanwidth = 1, timewidth = 1, stokes = 'I',
>>>          write='', MJDrange=[-1.0,-1.0], model=['delta'],
>>>          var=['p[0],p[1],p[2]'], p_ini=[0.0,0.0,1.0], phase_center = '',
>>>          fixed=[], fixedvar=[], scalefix='1.0', outfile = 'modelfit.dat',
>>>          NCPU=4, pbeam=False, ldfac = 1.22, dish_diameter=0.0, gridpix=0,
>>>          OneFitPerChannel=True, bounds=None, cov_return=False,
>>>          finetune=False, uvtaper=0.0, method='levenberg', wgt_power=1.0,
>>>          LMtune=[1.e-3,10.,1.e-5,200,1.e-3], SMPtune=[1.e-4,1.e-1,200],
>>>          only_flux=False, proper_motion = 0.0, HankelOrder = 80,
>>>          amp_gains = {}, phase_gains = {})

Let us assume that you have imported this module with the alias uvm (which is the case if you followed the installation instructions). To fit the flux density of a point source (centered at the coordinates of the field called M81), using the continuum emission in all the spws of a measurement set called myms.ms, just do: 

>>> myfit = uvm.uvmultifit(vis="myms.ms", spw="", model="delta", var="0,0,p[0]",
>>>         field="M81", OneFitPerChannel=False, outfile="results.dat")

The program will write the fitting result (and other useful metadata) into the file called results.dat. But it will also return a so-called uvmultifit instance (which, in this example, has been called myfit), that can be used to do a lot of things with your modelling. For instance, the fitting result (together with uncertainties, reduced Chi Square and covariance matrix, if the fit would have more than one fitting parameter) is stored in a dictionary that can be accessed as: 

>>> myfit.result

This dictionary has several keys worth noticing:

  • myfit.result['Parameters'] The fitting parameters. If the fit is done in spectral-line mode, these are organized per spw and channel.

  • myfit.result['Frequencies'] The frequency corresponding to each set of fitting parameters (useful for cases of fits in spectral-line mode).

  • myfit.result['Uncertainties'] The uncertainties of the fitted parameters.

    Note

    The uncertainties are estimated from the Jacobian matrix, and scaled so that the reduced Chi squared equals unity. Null (or ridicously small) uncertainties may indicate an unsuccessful fit. It is always a good idea to take a look at the post-fit residuals to assess the quality of your fit.

    Note

    The Jacobian is computed using numerical approximations of the derivatives of the fitting functions w.r.t. the parameters.

  • myfit.result['Reduced Chi Square'] The reduced Chi Squared, before rescaling the uncertainties.

    Note

    Notice that this is the reduced post-fit Chi squared, computed using the original visibility uncertainties (i.e., those read from the data). In other words, this is the reduced Chi Square computed before UVMultiFit scales the uncertainties to make the reduced Chi squared equal to unity. The user can check with this value whether the visibility uncertainties are well scaled to the natural scatter of the data.

    Note

    Quite large values of the reduced Chi Squared are likely indicative of too high data weights, which should be harmless to the fit, as long as the relative weights among visibilities are OK (in other words, a fit with UVMultiFit is insensitive to any global scaling factor of the weights).

    Note

    Too high values of the reduced Chi Squared could also be indicative of a bad fit. In such a case, though, the a-posteriori parameter uncertainties would also be high. The user should check the parameter uncertainties as an extra assessment of the quality of the fit.

  • myfit.result['covariance'] The full covariance matrix of the parameters (in case the user asked for it).

Note

There are other elements in your uvmultifit instance that may be interesting to you, if you want to do more advanced stuff. Look at the section Some Useful Methods for details.

If you made a fit in spectral-line mode and want to plot the first fitted parameter (i.e., p[0], see below for the syntax details) against frequency, for the third spectral window in the measurement set, the command could be (se assume that pylab has been imported): 

>>> pylab.plot(myfit.result['Frequency'][3], myfit.result['Parameters'][3][:,0])

To plot the second fitted parameter (i.e., p[1]), just execute: 

>>> pylab.plot(myfit.result['Frequency'][3], myfit.result['Parameters'][3][:,1])

and so on. If there is only 1 fitting parameter: 

>>> pylab.plot(myfit.result['Frequency'][3], myfit.result['Parameters'][3])

Note

The fitted parameters are ordered by spectral window in the order given by the “spw” parameter. Hence, if spw=‘2,3’, then the first element of the list of parameters (i.e., myfit.result['Parameters'][0]) will be the parameters fitted for the spw number 2.

Note

For fits to the continuum (i.e., all channels together, see below), myfit.result['Parameters'] will only have one entry (i.e., the parameter values fitted to all channels at once)

Model Examples

In the following examples, we only give the most important keywords to setup the fit. The rest of keywords (e.g., the name of the measurement set, the fields, the range of spectral windows, etc.) are not important to understand the model setup and are thus avoided for clarity.

Details about the variables of the models and how to set fitting parameters are given in the sections below.



  • EXAMPLE 0: A point source with free position and free flux density: 

    >>> model = ['delta']
>>> var = ['p[0],p[1],p[2]']

    In this case, the RA shift w.r.t. the image center is p[0] (in arcsec), the Dec shift is p[1] (also in arcsec), and the flux density is p[2] (in Jy). If we know that the source has to be located within 1 arcsec of the phase center, we can add this information as a set of bounds, i.e.: 

    >>> bounds = [[-1,1], [-1,1], None]

    If we also want to force the flux density to be positive: 

    >>> bounds = [[-1,1], [-1,1], [0,None]]


  • EXAMPLE 1: Two deltas, being the position of the second one fixed w.r.t. the position of the first one. Let’s say that the second delta is shifted at 0.5 and 0.6 arcsec, in RA and Dec (respectively), from the first delta, and we want the position of the first delta to be free in our fit. Then: 

    >>> model = ['delta', 'delta']
>>> var = ['p[0],p[1],p[2]', 'p[0]+0.5, p[1]+0.6, p[3]']

    In this case, p[0] and p[1] are the RA and Dec of the first delta; p[2] is the flux density of the first delta; and p[3] is the flux density of the second delta. Notice that the RA and Dec position of the second delta is equal to that of the first delta plus a fixed shift.



  • EXAMPLE 2: A ring plus a delta at its center. The absolute position of the compound source is free, and the ring is circular. 

    >>> model = ['ring','delta']
>>> var = ['p[0],p[1],p[2],p[3],1.0,0.0', 'p[0],p[1],p[4]']

    In this case, p[0] and p[1] are the RA and Dec of both components (i.e., the ring and the delta); p[2] is the total flux density of the ring and p[3] is its diameter; p[4] is the flux density of the delta. Notice that the axes Ratio of the ring is set constant (and unity) and the position angle is also set constant (although it’s meaningless for Ratio=1). For extended models, it is a good idea to bound the size to have positive values. Hence, in this case: 

    >>> bounds = [None, None, None, [0,None]]

    In case we also want to bound the fitted fluxes to be positive, we would have: 

    >>> bounds = [None, None, [0,None], [0,None], [0,None]]


  • EXAMPLE 3: Like Example 1, but fixing also the flux-density of the second delta to be 2.5 times that of the first delta: 

    >>> model = ['delta','delta']
>>> var = ['p[0],p[1],p[2]', 'p[0]+0.5, p[3]+0.6, p[2]*2.5']


  • EXAMPLE 4: A circularly-symmetric disc with a hole (i.e., with its inner half subtracted): 

    >>> model = ['disc','disc']
>>> var = ['p[0],p[1],4/3*p[2],p[3],1,0',
>>>        'p[0],p[1],-p[2]/3,p[3]/2,1,0']

    In this case, p[0] and p[1] are the RA and Dec shifts, respectively; p[2] is the flux density of the disc with the hole (i.e., with its inner half subtracted), and p[3] is the disc diameter. Notice that the hole in the disc has been created by adding a negative disc of size equals to 1/2 of the size of the larger disc, and flux density equals to -1/4 of that of the larger disc. The overall effect of both discs is that of one single disc with a hole (i.e., with no emission at radii < 1/2 of the outer radius).



  • EXAMPLE 5: A delta component with a spectral index, fitted to the whole dataset (i.e., setting OneFitPerChannel=False): 

    >>> model = ['delta']
>>> var = ['p[0],p[1],p[2]*(nu/1.e9)**p[3]']

    In this case, the fitted spectral index will be p[3], and p[2] will be the flux density at 1 GHz. Notice that p_ini (the list of initial values for the parameters) must also have the a priori value of the spectral index. For this example, p_ini could be (for a source close to the center, with an approximate flux of 2.3 Jy at 1 GHz and an a priori spectral index of -0.7): 

    >>> p_ini = [0.0, 0.0, 2.3, -0.7]

    If the spectral index is well known, it can of course be fixed in the fit. In such a case: 

    >>> model = ['delta']
>>> var = ['p[0],p[1],p[2]*(nu/1.e9)**(-0.7)']
>>> p_ini = [0.0,0.0,2.3]

    Note

    Fitting sizes and spectral indices at the same time may produce crazy results, since both quantities are quite coupled in Fourier space. Hence, some care should be taken when fitting to all frequencies together. Notice also that, unless your bandwidth is quite wide and/or your SNR is quite high, any fit to the spectral index may not be reliable.



  • EXAMPLE 6: A filled sphere with its inner half (i.e., the core) removed. The sphere is fixed at the image center: 

    >>> model = ['sphere','sphere']
>>> var = ['0,0,9/8*p[0],p[1],1,0','0,0,-p[0]/8,p[1]/2,1,0']

    In this case, p[0] is the total flux density and p[1] is the outer diameter. This example is similar to that of the disc with a hole, but in this case the model is a sphere, where we have removed all the emission at radii < 0.5 times the outer radius.



  • EXAMPLE 7: A disc with a hole of variable size. In this case, the smaller disc with negative flux density (which is used to generate the hole) must always have the same surface brightness (in absolute value) that the larger positive disc. Hence, and if we fix the disc position at the image center (for simplicity), we have: 

    >>> model = ['disc','disc']
>>> var = ['0,0,p[0]*p[1]**2./(p[1]**2.-1),p[2],1,0',
>>>        '0,0,-p[0]/(p[1]**2.-1),p[2]/p[1],1,0']

    In this case, the flux density of the disc with the hole is p[0], the outer dimaeter is p[2], and p[1] is the ratio of the outer size to the inner size.



  • EXAMPLE 8: Simultaneous self-calibration and source fitting. We fit a source flux density and perform Global Fringe Fitting (i.e., we fit for phases, delays, and delay rates) in one shot. Let’s fit stokes RR and assume that there are 3 antennas (for simplicity): 

    >>> stokes = 'RR'
>>> model = ['delta']
>>> phase_gains = {1:'p[0] + 6.2832*(p[1]*(nu-nu0) + p[2]*t)',
>>>                2:'p[3] + 6.2832*(p[4]*(nu-nu0) + p[5]*t)'}
>>> var = '0,0,p[6]'

    Good initial estimates for the phase gains can be obtained using the “QuinnFF” method.

Model Details

Each model depends on a list of variables that can be written as any algebraic combination of fitting parameters. Some explanatory examples are given in the section Model Examples (see also help(uvm), if you imported the module with the uvm alias).

Here is the list of currently-supported models, and their variables:

  • delta -> Variables: RA, Dec, Flux
  • disc -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • ring -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • Gaussian -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • sphere -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • bubble -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • expo -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • power-2 -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • power-3 -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle
  • GaussianRing -> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle, Sigma

These models have the following meaning:

  • sphere stands for an optically-thin uniform filled sphere.

  • bubble stands for a uniform spherical surface.

  • expo stands for an exponential radial flux decay.

  • power-2 stands for a decay as \(1/(r^2 + r_0^2)\) (notice that in this case, the flux is the integral from r=0 to r=r0)

  • power-3 stands for a decay as \(1/(1 + (2^{2/3} - 1)(r/r_0)^2)^{3/2}\).

  • GaussianRing stands for a radial profile following the expression: \(F*exp(-(Q - R)^2/(2*Sigma^2))\)

    Note

    The GaussianRing is built by a series expansion of the Bessel function in the computation of the Hankel transform. The order of the expansion (default is 80) can be set in the new HankelOrder keyword. The larger the source (in beam units) the higher HankelOrder should be, to keep the numerical precision.

The meaning of the model variables is:

  • RA and Dec are the shifts w.r.t. the phase center (in arcsec)

  • Flux is the total flux density of the component (in Jy)

  • Major is the diameter along the major axis

  • Ratio is the size ratio between the reference axis and the other axes (i.e., it is set to 1.0 for circularly-symmetric sources).

    Warning

    If Ratio is forced to be higher than 1.0, then Major (see above) will indeed become the minor axis! If you want to avoid this issue, you should bound Ratio to be lower than (or equal to) 1.0. Otherwise, you will have to remember this issue when you interprete your fitted parameters.

  • PositionAngle is the angle of the reference axis, from North to East (in deg.)

  • Sigma (whenever it is used) is an auxiliary variable for models that need more than one ‘size-like’ parameter to be defined (e.g., the ‘GaussianRing’, where we have the size of the ring and its width).

UVMultiFit Data Properties

UVMultiFit will return an instance with many properties and methods. Some examples are listed below.

  • averdata: A set of arrays (one per spectral window) with the visibilities, averaged in time and/or frequency (according to the chanwidth and timewidth parameters).
  • averweights: The (square root of the) weights used in the fit.
  • u, v, and w: the coordinates in Fourier space of the visibilities (one array per spectral window).
  • t: the observing times (Modified Julian Dates). One array per spw.

Note

A UVMultiFit instance can occupy quite a bit of physical memory. It may be a good idea to delete the instance from memory (i.e., using del myfit and gc.collect()) or restart CASA once the user has got the desired results from the fit (this is the recommended approach).

The most important methods of UVMultiFit are described below.

The UVMultiFit Class

class uvmultifit.uvmultifit(vis='', spw='0', column='data', field=0, scans=[], uniform=False, chanwidth=1, timewidth=1, stokes='I', write='', MJDrange=[-1.0, -1.0], ldfac=1.22, model=['delta'], var=['p[0],p[1],p[2]'], p_ini=[0.0, 0.0, 1.0], phase_center='', fixed=[], fixedvar=[], scalefix='1.0', outfile='modelfit.dat', NCPU=4, pbeam=False, dish_diameter=0.0, gridpix=0, OneFitPerChannel=True, bounds=None, cov_return=False, finetune=False, uvtaper=0.0, method='levenberg', wgt_power=1.0, LMtune=[0.001, 10.0, 1e-05, 200, 0.001], SMPtune=[0.0001, 0.1, 200], only_flux=False, proper_motion=0.0, HankelOrder=80, phase_gains={}, amp_gains={})

This is the main UVMultiFit class.

It includes pointers to the data and many different methods to perform fits, retrieve Chi square values, covariance matrices, calibrate data, write either model or residual visibilities to your measurement sets, etc.

Fits with different models can be performed on the same data, (or on different time ranges within the same data to, e.g., perform variability studies) with no need to re-read the same data every time.

The multi-threaded C++ engine, coupled to the flexibility of Python and CASA, makes UVMultiFit a fast, flexible and powerful tool for your advanced visibility-based analysis.

Parameters:
vis : str
Name of the measurement set. It can also be a list of MS names.
spw : str
String (or list of strings). These are the spectral window(s) and channel selection(s) to be fitted. For each spw, the user can select many channel ranges in the usual CASA way. If one string is given, all ms listed in ‘vis’ are supposed to have the same spectral configuration. If a list of strings is given, one per ms, this restriction doesn’t apply.
column : str
The data column. It can be either ‘data’ or ‘corrected’
field : str
The id number (or name) of the target source in the measurement set(s).
pbeam : bool
If false, the primary-beam correction is not applied. This is very important for fitting mosaic data.
dish_diameter : double
In case that the antenna diameters cannot be read from the datasets, the user must provide the antenna diameters (in meters).
ldfac : double
Proportionality constant between the ratio ‘lambda/Diameter’ and the FWHM of the primary beam (assumed to be a Gaussian!). I.e.: FWHM = ldfac*lambda/Diameter. Normally, ldfac = 1.22 should be fine, although 1.00 works better with data coming from simobserve.
scans : list
List of integers; default []. The id numbers of the scans to load. Default means to load all the scans of the selected source. If multiple measurement sets are selected, this should be a list of lists (i.e., one list of scans per measurement set). For instance, if vis=["ms1.ms","ms2.ms"], then scans=[[1,2,3],[]] would select the scans 1, 2, and 3 from the ms1.ms dataset and all the scans from the ms2.ms dataset.
chanwidth : int
Number of spectral channels to average in each chunk of data before the fit. BEWARE of the bandwidth-smearing effects if your baselines are long (and/or your source is large).
timewidth : int
Number of time channels (i.e., integration times) to average in each chunk of data before the fit. The averaging is always cut at the scan boundaries (so a very large timewidth means to average down to one visibility per scan). BEWARE of the time-smearing effects!
MJDrange : list
List of two floats. These are the initial and final Modified Julian Dates of the data to be used in the fitting. Notice that all the scans asked by the user are loaded a-priori, and the MJDrange condition is applied afterwards. This way, variability studies can be performed efficiently, by loading all data at once and then setting different MJDranges iteratively. Default (i.e., <= 0.0) means not to select data based on JD time range.
stokes : str
Polarization product. Can be any of PI, I, Q, U, V. It also accepts individual correlation products: XX, YY, XY, YX, RR, LL, LR, or LR. Default is I. If PI (which stands for Polarization Independent) is given, the program will compute I whenever possible and use either XX, YY, RR or LL otherwise. This way, the amount of data used in polarization-independent fits is maximized.
model : list
List of strings (i.e., model components to fit). Each component is given as a string. Possible models are: delta, disc, Gaussian, ring, sphere, bubble, expo, power-2, power-3, and GaussianRing. If only one model component is being fitted, the model keyword can also be a string.
var : list
List of strings (or just one string, if only one model component is being fitted). These are the variables for each model. The variables can be set to any algebraic expression involving the fitting parameters (being the ith parameter represented by p[i]) and the observing frequency in Hz (represented by nu). Any numpy function can also be called, using the prefix np (e.g., p[0]*np.power(nu,p[1])). See some examples below.
p_ini : list
List of the initial values of the fitting parameters. This is expected to be a list of floats.
phase_center : str
The sky position where all components are referenced to. If an empty string is given, the phase center will be that of the first field id that is being read (i.e., if a mosaic is being read, the first pointing will be set as the phase center). If the string is not empty, the program expects to find a coordinate in CASA format (i.e., J2000 RA Dec, where RA is in format 00h00m00.0s and Dec is in format 00d00m00.0s).
fixed : list
Like model, but defines model components with completely fixed variables (i.e., whose variables are defined only by numbers; not fitting parameters). This model will be computed only once (i.e., just before the fit), hence making the code execution faster. The user can load the model column of the measurement set(s) as a fixed model, by setting fixed='model_column'.
fixedvar : list
Like var, but refers to the fixed model. Hence,it is expected to be either a list of numbers or a list of strings representing numbers. This is not needed if fixed = 'model_column' (since, in that case, the model column is read as is from the measurement set).
scalefix : str
String representing a function that sets a scaling factor for the fixed-model’s total flux density. It can be a function of the fitting parameters (e.g., scalefix='p[0]' will multiply the overall flux density of the fixed model by p[0]) and can also be a function of the observing frequency nu.
OneFitPerChannel : bool
If True, independent fits are performed to the different frequency channels, one by one. If False, one common fit is performed to all data. In this case, the user may want to fit for the spectral indices of the components, if the fractional bandwidth is wide.
outfile : str
Name of the output file to store results (i.e., fitting parameters, uncertainties, and metadata, in ascii format).
bounds : list
List of boundaries (i.e., minimum and maximum allowed values) for the fitting parameters. ‘None’ means that no bound is defined. If the list is empty, no bounds are assumed for any parameter (see examples below).
cov_return : bool
If True, the covariance matrix for each fit is added to the dictionary returning from the fit() method (the dictionary key will have the name 'covariance', see the Returns section below).
uvtaper : double
Default is 0.0. If not 0.0, the weights of the visibilities are multiplied by a Gaussian in Fourier space, whose HWHM is the value of uvtaper, in meters.
uniform : bool
Default is False. If True, the weights of all data are made equal. Notice that a uvtaper can still be applied (i.e., by setting the uvtaper parameter).
finetune : bool
Default is False. If set to True, the fit is not performed, but only a uvmultifit instance is created with the data properly read and the models compiled and ready. The user can then run different methods of the UVMultiFit class by him/herself (see the help text for each method) before actually fitting the data. This can be useful, for instanse, if the user wants to try many different models (and/or subtract many different fixed models), apply ad-hoc calibrations, perform an MCMC exploration of the parameter space, etc., without having to reload the data every time after each step.
wgt_power : double
Default is 1. Power index of the visibility weights in the computation of the Chi square. wgt_power = 1 would be the statistically justified value, although other values may be tested if the user suspects that the visibility uncertainties are not well estimated in his/her dataset.
method : str

Method to use in the chi-square minimization. Default is simplex. Possible values are simplex and levenberg. Sometimes, the least-squares minimization may not converge well with the Levenberg-Marquardt method (if the model is complicated and/or the uv coverage is quite sparse). Levenberg-Marquardt also requires a lot of memory, so may not be very convenient if the datasets are very large. In these cases, a Chi square minimization using the simplex algorithm may work beter. However, simplex may require more function evaluations to find the minimum.

Warning

The SIMPLEX method does not currently provide parameter uncertainties. It just returns the reduced Chi squared!

write : str
The kind of information to store in the measurement set(s) after the fit. Default is ‘’ (i.e., does not change anything in the datasets).

  • If it is set to model, the best-fit model is saved in the model column of the measurement set.

  • If it is set to residuals, the fit residuals are saved in the corrected’ column of the measurement set.

  • If it is set to calibrated (this option will be available in the next release), the gains defined in the amp_gains and phase_gains dictionaries (see below) will be applied, and the calibrated data will be saved in the corrected column of the ms.

    Currently, this keyword is only activated if stokes is set to either PI, I or an individual correlation product (like XX or XY) and if both timewidth and chanwidth are set to 1.

NCPU : int
Default is 4. Number of threads allowed to run in parallel.
SMPtune : list
Used to fine-tune the Simplex algorithm. Change only if you really know what you are doing. The meaning of the list elements is:
  • SMPtune[0] -> Maximum allowed error in the parameters, from the search of the Chi Square minimum. Default is 1.e-4.
  • SMPtune[1] -> Relative size of the first step in the parameter search. Default is 1.e-1.
  • SMPtune[2] -> Maximum number of iterations allowed per fitting parameter. Default is 200
LMtune : list
Used to fine-tune the Levenberg-Marquardt algorithm. Change only if you really know what you are doing. The meaning of the list elements is:
  • LMtune[0] -> The Lambda factor to weight up the Newton component (i.e., H + Lambda*H_diag = Res), where H is the Hessian and Res the vector of residuals. Default is 1.e-3
  • LMtune[1] -> The factor to multiply (divide) Lambda if the iterator worsened (improved) the Chi Squared. Default: 10.
  • LMtune[2] -> The maximum relative error allowed for the ChiSq. Default: 1.e-5
  • LMtune[3] -> Maximum number of iterations allowed per fitting parameter. Default: 200
  • LMtune[4] -> (optional) Maximum relative error allowed for the parameters. Default: 1.e-3. If it is not provided, then LMtune[4] = LMtune[2]
only_flux : bool
Default is False. If True, the program assumes that only the flux densities of all components are going to be fitted. Furthermore, the fitting parameters shall be just equal to the flux densities being fitted, and should be given in the same order as the list of model components. In these cases, using only_flux = True can speed up the fit quite a bit, especially if there are many components to be fitted. This option may be useful to implement simple compressed sensing-like algorithms.
proper_motion : list
List of 2-element lists of numbers: Each element (i.e., each list of two numbers) is the proper motion, in RA and Dec, of each model component. The units are arc-seconds per year. Proper motions cannot be fitted yet, but may be fittable in future versions of UVMultiFit. Default is a float = 0.0, meaning that all proper motions are null. If the proper motions are not null, the position used as reference for the fit will correspond to that of the first integration time of the first scan of the observed field.
HankelOrder : int
Only used for models without an analytic expression (i.e., at the moment, only the GaussianRing model). In these cases, UVMultiFit performs the Hankel transform by using the series expansion of the Bessel function J0. The order of this expansion is set by this keyword. The larger the distance in Fourier space (i.e., the larger the source, in beam units), the higher HankelOrder should be, to keep the numerical precision. Default is 80. For cases of sources with sizes similar to the synthesized beam, HankelOrder=80 should suffice. The use of too high values may cause Overflow errors!
amp_gains: dict
Dictionary (default empty). Controls whether to solve for antenna amplitude gains and source-model parameters simultaneously. The keys of the dictionary are the indices of the antennas to self-calibrate. The values of the dictionary are strings, whose elements represent functions of the fitting parameters. See examples below.
phase_gains : dict
Same as for amp_gains, but dealing with phase gains.

Some Useful Methods

uvmultifit.checkInputs()

Reads all the inputs parameters and performs sanity checks.

You have to run this everytime after changing any parameter of the UVMUltiFit instance.

Many self-consistency checks are performed, related to data selection, number of internal variables reading, averaging, and formatting, etc. It is always a good idea to run this method before fitting, if the user has changed some parameters (to redo a fit) or is using finetune=True.

Warning

If you change any keyword related to the data selection (with the exception of MJDRange) you must run readData() after checkInputs(), to ensure that the new data are properly read and set before the fit.

The keywords that need a rerun of readData() are: vis, spw, stokes, column, scan, field, uniform, timewidth, chanwidth, phase_center, wgt_power, uvtaper.

Warning

If you change the equations of your fitting model (ond/or those of the antenna gains) you must run initModel() before the fit.

The keywords that need a rerun of initModel() are: model, var, fixed, fixedvar, scalefix

Notice that just changing p_ini or bounds is OK. You do not need to reinit the model to repeat the fit.

uvmultifit.readData(del_data=True)

Reads the data, according to the properties vis, column, chanwidth, etc.

It then fills in the properties averdata, averfreqs, averweights, u, v, w, etc.

Each one of these properties is a list with the data (one list item per spectral window/scan).

A previous successful run of function checkInputs() is assumed.

Note

Instead of re-reading data from scratch, using the same uvmultifit instance, a better approach may be to restart CASA and create a fresh uvmultifit instance with the new data, avoiding some memory leakage related to potential hidden references to the data in the IPython’s recall prompt.

uvmultifit.initModel()

Allocates memory for the modeler data, which will be used by the C++ extension.

Also compiles the model variables. It is a good idea to run this method everytime that the user reads new data to be fitted (i.e., after a call to readData())

It is indeed MANDATORY to run this function if the model to be fitted has changed.

uvmultifit.fit(redo_fixed=True, reinit_model=False, save_file=True, nuidx=-1, reset_flags=False)

Fits the data, using the models previously compiled with initModel().

reinit_model : bool
It is False by default. If False, the models used are those already compiled. If True, the models are recompiled, according to the contents of the model, var, fixed, fixedvar properties, and all the references to the data arrays are refreshed.
redo_fixed : bool

It is True by default. If False, the fixed model will not be recomputed throughout the fit.

Warning

Setting redo_fixed=False can be dangerous if you are fitting in spectral-line mode and have channels with very different frequencies (i.e., a wide fractional bandwidth), since the UV coordinates will not be re-projected in that case.

save_file : bool
It is True by default. If False, the external ascii file with the fitting results will not be created.
nuidx : int
A helper parameter for internal use (if -1, fit to the continuum. Otherwise, fit only the nui channel). The user should NOT need to redefine (or force) this parameter when calling this function.
reset_flags : bool
Default is False. This is used to clear out the status of bad data that may have been set by spetial routines of UVMultiFit (e.g., the Quinn Fringe Fitter). Default is to NOT reset flags (this option should be OK most of the time).
uvmultifit.writeModel()

Writes the requested information into the measurement sets.

The information can be either the predictions of the compiled model(s) (i.e., they are written into the model column of the measurement set(s)), or the post-fit residuals (they are written into the corrected column) or the calibrated data (they are written into the corrected column as well). The actual information to write is set by the value of the write keyword of the UVMultiFit instance when the fit() method was called.

This function is executed only if the stokes keyword is set to either PI, I or an individual correlation product (like, e.g., XX or XY) and if no averaging has been performed neither in time nor frequency (i.e., if both timewidth and chanwidth are set to 1). This function should be called AFTER having run fit().

uvmultifit.chiSquare(p)

Returns a list with 2 items: The Chi Square value, computed for the parameter values given in the p list) and the number of degrees of freedom.

Parameters:
p : list
List of parameter values where the Chi square will be computed.
uvmultifit.QuinnFF(IF, refant, doModel, doGlobal)

Perform a Global Fringe Fitting search in delay-rate space.

It uses the Quinn estimator for the fringe peak.

Parameters:
IF : int
The spectral window to fringe-fit (index referred to the data already read).
refant : int
The index of the antenna to use as reference. If this antenna is missing or has bad data, another refant will be picked automatically.
doModel : bool
If True, deconvolve the fixed model before the GFF.
doGlobal : bool
If True, globalize the solutions (i.e., use all available baselines to estimate the antenna gains).
Returns:

Returns a list of 5 elements:

  • An array with the delays (in seconds). One value per antenna (the reference antenna, as well as any other antenna with failed solutions, will have null values).
  • An array with the rates (in 1/s), one value per antenna.
  • An array with the phases (in radians), one value per antenna.
  • A double (the delay equivalent to one pixel in the delay/rate matrices).
  • A double (the rate equivalent to one pixel in the delay/rate matrices).

The IMMultiFit Class

class uvmultifit.immultifit(parent=None, reinvert=True, start=0, nchan=-1, psf='', residual='', dBcut=-30.0, **kwargs)

Similar to uvmultifit, but the engine uses images created by the CASA task clean.

All keywords are equal to those of uvmultifit, with the exception of these:

Parameters:
psf : str
The image (or image cube) of the Point Spread Function. Usually, the name of this image is in the form blablabla.psf
residual : str

The image (or image cube) of the clean residuals. Usually, the name is in the form blablabla.residual.

Note

You should clean with niter=0, to ensure that no signal is deconvolved.

Note

If you want primary-beam correction, you should use the blablabla.image, instead of blablabla.residual.

reinvert : bool
If immultifit has already been run and the user wants to refit the same data, the Fourier inversion of the images (or image cubes) can be avoided if reinvert=False. Default is True
dBcut : double
Mask the uv points with a signal lower than dBCut of the peak (in dB). Use this to minimize convolution-like artifacts in UV space. Default: -30dB => SNR = 1000.

Note

The keywords vis, spw, timewidth, chanwidth, phase/amp_gain, etc. are not used by immultifit.

Note

Instead of writing models or residuals into a measurement set, immultifit writes images with the same gridding as the original images used.