UVMultiFit’s Documentation¶
A CASAbased flexible visibilityfitting 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 multicomponent models to the visibilities in one (or several) measurement set(s). The fits can be performed to the continuum or in spectralline mode. Advanced model structures can be implemented easily and several different fits can be performed with no need of reloading the data.
Installation¶
Just follow these steps:
 Download the latest version from Launchpad.

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.
 Install the dependencies (i.e., GNU Scientific Library, so far).
In GNU/Linux, execute one of these commands:
sudo aptget install libgsl2 sudo aptget install libgsldev sudo aptget install libgsl0ldblIn MacOS, execute all these commands:
sudo port install gsl export LIBRARY_PATH="/opt/local/lib" export LD_LIBRARY_PATH="/opt/local/lib"

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

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 executecd ~ ; pwd
in a terminal).

NoteIf 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==CASAVERSIONOFNUMPY).
 Enjoy!
Any feedback, problem with the installation/running and/or bug report should be sent either to the ARC Nordic Node (contact@nordicalma.se) or to the source maintainer (ivan.martividal@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.e3,10.,1.e5,200,1.e3], SMPtune=[1.e4,1.e1,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 socalled 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 spectralline 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 spectralline 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 postfit 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 postfit 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 aposteriori 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 spectralline 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 isp[1]
(also in arcsec), and the flux density isp[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]
andp[1]
are the RA and Dec of the first delta; p[2] is the flux density of the first delta; andp[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]
andp[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 andp[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 forRatio=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 fluxdensity 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 circularlysymmetric 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]
andp[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), andp[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]
, andp[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 andp[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 isp[2]
, andp[1]
is the ratio of the outer size to the inner size.

EXAMPLE 8: Simultaneous selfcalibration 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]*(nunu0) + p[2]*t)', >>> 2:'p[3] + 6.2832*(p[4]*(nunu0) + 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 currentlysupported models, and their variables:
delta
> Variables: RA, Dec, Fluxdisc
> Variables: RA, Dec, Flux, Major, Ratio, PositionAnglering
> Variables: RA, Dec, Flux, Major, Ratio, PositionAngleGaussian
> Variables: RA, Dec, Flux, Major, Ratio, PositionAnglesphere
> Variables: RA, Dec, Flux, Major, Ratio, PositionAnglebubble
> Variables: RA, Dec, Flux, Major, Ratio, PositionAngleexpo
> Variables: RA, Dec, Flux, Major, Ratio, PositionAnglepower2
> Variables: RA, Dec, Flux, Major, Ratio, PositionAnglepower3
> Variables: RA, Dec, Flux, Major, Ratio, PositionAngleGaussianRing
> Variables: RA, Dec, Flux, Major, Ratio, PositionAngle, Sigma
These models have the following meaning:

sphere
stands for an opticallythin uniform filled sphere. 
bubble
stands for a uniform spherical surface. 
expo
stands for an exponential radial flux decay. 
power2
stands for a decay as \(1/(r^2 + r_0^2)\) (notice that in this case, the flux is the integral fromr=0
tor=r0
) 
power3
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 circularlysymmetric 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 ‘sizelike’ 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 thechanwidth
andtimewidth
parameters).averweights
: The (square root of the) weights used in the fit.u
,v
, andw
: 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, 1e05, 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 reread the same data every time.
The multithreaded C++ engine, coupled to the flexibility of Python and CASA, makes UVMultiFit a fast, flexible and powerful tool for your advanced visibilitybased 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 primarybeam 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"]
, thenscans=[[1,2,3],[]]
would select the scans 1, 2, and 3 from thems1.ms
dataset and all the scans from thems2.ms
dataset.  chanwidth : int
 Number of spectral channels to average in each chunk of data before the fit. BEWARE of the bandwidthsmearing 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 timesmearing 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 apriori, 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
, orLR
. Default isI
. IfPI
(which stands for Polarization Independent) is given, the program will computeI
whenever possible and use eitherXX, YY, RR
orLL
otherwise. This way, the amount of data used in polarizationindependent 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, power2, power3
, andGaussianRing
. 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 bynu
). Any numpy function can also be called, using the prefixnp
(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 format00h00m00.0s
and Dec is in format00d00m00.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 fixedmodel’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 byp[0]
) and can also be a function of the observing frequencynu
.  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 adhoc 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 chisquare minimization. Default is
simplex
. Possible values aresimplex
andlevenberg
. Sometimes, the leastsquares minimization may not converge well with the LevenbergMarquardt method (if the model is complicated and/or the uv coverage is quite sparse). LevenbergMarquardt 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 bestfit 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 (likeXX
orXY
) 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 finetune 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.e4.SMPtune[1]
> Relative size of the first step in the parameter search. Default is 1.e1.SMPtune[2]
> Maximum number of iterations allowed per fitting parameter. Default is 200
 LMtune : list
 Used to finetune the LevenbergMarquardt 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.e3LMtune[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.e5LMtune[3]
> Maximum number of iterations allowed per fitting parameter. Default: 200LMtune[4]
> (optional) Maximum relative error allowed for the parameters. Default: 1.e3. If it is not provided, thenLMtune[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 sensinglike algorithms.  proper_motion : list
 List of 2element 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 arcseconds 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 sourcemodel parameters simultaneously. The keys of the dictionary are the indices of the antennas to selfcalibrate. 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 selfconsistency 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 runreadData()
aftercheckInputs()
, 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
orbounds
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 rereading 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 spectralline mode and have channels with very different frequencies (i.e., a wide fractional bandwidth), since the UV coordinates will not be reprojected 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 postfit 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 thefit()
method was called.This function is executed only if the
stokes
keyword is set to eitherPI
,I
or an individual correlation product (like, e.g.,XX
orXY
) and if no averaging has been performed neither in time nor frequency (i.e., if bothtimewidth
andchanwidth
are set to 1). This function should be called AFTER having runfit()
.

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 delayrate space.
It uses the Quinn estimator for the fringe peak.
Parameters:  IF : int
 The spectral window to fringefit (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 primarybeam correction, you should use the
blablabla.image
, instead ofblablabla.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 convolutionlike 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.