Fit an arbitrary sky (image) model to an existing MS

Fitting a Sky Model to an Existing MS in the UV Domain

20may2015 CASA v4.3.1

This casaguide illustrates how to fit an arbitrary sky image model to an existing MS. Along the way it illustrates how to "get your hands" on the raw visibility data in the python interface and how to compute the Chi-squared of the visibility data to your given model. This could be useful if you have a set of source models you would like to compare in the UV domain, but the models are too complicated to be handled by the set of models provided in CASA's uvmodelfit() task. A couple of ALMA 7-m specific points are illustrated along the way.

# import FITS image to a CASA image
#  note: the FITS image need not have single
#    precision data unless defaultaxes=True (but we are
#    using the simultil.modifymodel() method to fix
#    the axies, so that's unnecessary)
#fitsimage='lyszmapjyperpixRightSignSingPrec.fits'
default importfits
fitsimage='lyszmapjyperpixRightSign.fits'
imagename='lyszmapjyperpixRightSign.image'
overwrite=True
defaultaxes=False
inp
go

from simutil import simutil
util=simutil()

#
# i think 2GHz for inwidth was an illegal spectral unit
#  ? so changed to MHz.
# ****NOTE - the incenter will have to be chosen to be consistent
#      for each SPW***
util.modifymodel(inimage='lyszmapjyperpixRightSign.image',
                 outimage='lyszmapjyperpixRightSignFixed.image',
                 inbright='',indirection='',incell='',incenter='85.015425GHz',
                 inwidth='1937.5MHz',innchan=1)


# split out one SPW and the field(s) of interest-
#  make two copies, one of which will get overwritten by the MOCK
spwid='0'
realmsname='calibrated_final_cont_spw'+spwid+'_AVG.ms'
split(vis='calibrated_final_cont.ms',outputvis=msname,field='4',
      spw=spwid,combine='scan',timebin='1e8')
msname='calibrated_final_cont_spw'+spwid+'_AVG_MOCK.ms'
split(vis='calibrated_final_cont.ms',outputvis=msname,field='4',
      spw=spwid,combine='scan',timebin='1e8')
#old form -
#  split(vis='calibrated_final_cont_SAVE.ms',outputvis=msname,field='4',spw=spwid)

#
# set the telescope name manually- the need for this is probably a defect -
#  (sim tool infers primary beam from telescope name)
# ***DO NOT DO THIS FOR THE 12M DATA!!!***
#
tb.open(msname+'/OBSERVATION',nomodify=False)
tel=tb.getcol("TELESCOPE_NAME")
for i in range(len(tel)):
  tel[i]="ACA"

tb.putcol("TELESCOPE_NAME",tel)
tb.done()

sm.openfromms(msname)

#sm.setdata(spwid=0,fieldid=4)
#sm.setdata()
# this sets the primary beam
sm.setvp()
sm.summary()

sm.predict(imagename='lyszmapjyperpixRightSignFixed.image')
# ***NOTE: the result has gone in the data column, not the model column

# check result with sanity -
plotms(vis=msname,xaxis='uvwave',yaxis='amp',ydatacolumn='data',avgtime='1e8')

# plot what the observed model would look like, sans noise-
clean(vis=msname,niter=0,cell='2.5arcsec',imsize=[128,128],imagename='model7m')

#
# issues with sim toolkit
#  -need to do one SPW at a time. pre-split it.
#  -by default it seems to be using the 12m PB for the 7m data.
#    (and the 12m is a rough approx. to the actual 12m beam,
#    ie it is an airy disk rescaled to have ~ same main lobe size
#    as the real beam) --> fix by tweaking OBSERVATION table
#    for ACA it then uses 6.5m airy disk.
#

###############################

# extract visibilities to the python layer
#  using the measurement set tool
ms.open(msname)
# the selection can be tricky, but it seems to work fine with
#  the default parameters using the SPLIT() strategy we have
#  done above-
ms.selectinit()
# ms.getdata() returns a dictionary, each element
#  of which is a numpy array
mydat=ms.getdata(['real','imaginary','u','v','weight'])
uu=mydat['u']
vv=mydat['v']
rr=mydat['real']
ii=mydat['imaginary']
# i got the weight from the sim dataset,should get
#  from real dataset -- but really should be identical--
wt=mydat['weight']
ms.close()

# now get visibility values from the actual (not mock) data-
ms.open(realmsname)
ms.selectinit()
realdat=ms.getdata(['real','imaginary','u','v','weight'])
r2=realdat['real']
i2=realdat['imaginary']
ms.close()

# what are the dims of the data?
r2.shape

# built-in sqrt object does not work elementwise
#  so import numpy's version -
import numpy as np
uvrad=np.sqrt(uu**2+vv**2)
ampl=np.sqrt(rr**2+ii**2)
a2=np.sqrt(r2**2+i2**2)

# first index indexes polarization (there are two --XX, YY)
# second is null (probably SPW or IF, nominally)-
pl.plot(uvrad,ampl[0,0,:],'.')
# hey, that looks familiar!

# calculate chi2-
c2= (r2[:,0,:]-rr[:,0,:])**2 * wt + (i2[:,0,:]-ii[:,0,:])**2 * wt
# typical data point contributes 1.6 to chi2-
np.median(c2)
# total chi2 per DOF-
c2.sum()/c2.size
# = 2.31 - better than i expected!

# nulling the model gives delta(chi2) = 51.4
#  -that's delta total chi2 for 36,736 DOF

# simple single-parameter linear least squares fit of the model to the data-
numerator =  r2[:,0,:] * rr[:,0,:] * wt + i2[:,0,:]*ii[:,0,:] * wt
denominator =  (rr[:,0,:])**2 * wt + (ii[:,0,:])**2 * wt
# LLS scale factor between data and model -
scalefac= numerator.sum() / denominator.sum()
# uncertainty in scalefac, if the weights are correct
#  in an absolute sense -
escalefac=sqrt(1.0/denominator.sum())