################################################                                             #
# Regression/Benchmarking Script for NGC 1333 #
#                                             #
###############################################

import time
import os

os.system('rm -rf ngc1333*.ms n1333.ms src.task.image n1333*.*.*cal src.t*.ms gcal.t*.ms src.tmom* src.t*.resid src.t*.model')


restore()

startTime = time.time()
startProc = time.clock()

datapath=os.environ.get('AIPSPATH')+'/data/regression/ATST2/NGC1333/'

print '*** 02 MAY ***'
print '--Import--'
###################################################
## Fill data from UVFITS format to measurementset
default('importvla')
importuvfits(fitsfile=datapath+'N1333_1.UVFITS',vis='ngc1333_1.ms')
importtime = time.time() 
##################################################
#print '--Observation summary--'
#listobs('ngc1333_.ms')
#listtime = time.time()
##################################################
## Flag Autocorrelations
print '--Flag Autocorrelations--'
default('flagautocorr')
flagautocorr('ngc1333_1.ms')
##################################################
###Flag antenna with antennaid 8 
default('flagdata')
flagdata('ngc1333_1.ms',antennaid=8,fieldid=-1)
#################################################
###flag all data whose amplitude  are not in range [0.0,2.0] on the
###parallel hands
default('flagdata')
flagdata('ngc1333_1.ms',clipcorr='RR',clipminmax=[0.0,2.0],fieldid=-1)
default('flagdata')
flagdata('ngc1333_1.ms',clipcorr='LL',clipminmax=[0.0,2.0],fieldid=-1)
default('flagdata')
#flagdata('ngc1333_both.ms',antennaid=13,timerange=['02-MAY-2003/18:30:00','02-MAY-2003/19:30:00'],fieldid=-1)
#default('flagdata')
#flagdata('ngc1333_both.ms',chans=[59,60,61,62],fieldid=-1,spwid=-1)
####################################################
### Flag data (which is bad) in  time range 
default('flagdata')
flagdata('ngc1333_1.ms',timerange=['02-MAY-2003/21:50:40','02-MAY-2003/22:01:10'],fieldid=-1)
###################################################
### Flag all antenna 13 data in the time range stated
default('flagdata')
flagdata('ngc1333_1.ms',antennaid=13,timerange=['02-MAY-2003/18:50:50','02-MAY-2003/19:13:30'],fieldid=-1)

flagtime=time.time()

#####set the model flux density of calibrater
print '--Setjy--'
default('setjy')
setjy('ngc1333_1.ms',field='0542+498_1') 
setjy('ngc1333_1.ms',field='0542+498_2') 
setjytime = time.time()

########################################
#####derive the gain calibration solutions
print '--Gaincal--'
default('gaincal')
gaincal('ngc1333_1.ms','n1333.1.gcal',mode='channel',nchan=55,start=4,msselect='FIELD_ID IN [0,12,14] && DATA_DESC_ID IN [0]',opacity=False,solint=0.,refant='27')
gaincal('ngc1333_1.ms','n1333.2.gcal',mode='channel',nchan=55,start=4,msselect='FIELD_ID IN [6,13,15] && DATA_DESC_ID IN [1]',opacity=False,solint=0.,refant='27')
gaintime = time.time()

#####################################
###derive the bandpass solutions
print '--Bandpass--'
default('bandpass')
bandpass('ngc1333_1.ms','n1333.1.bcal',msselect='FIELD_ID==14 && DATA_DESC_ID==0',opacity=False,gaintable='n1333.1.gcal',refant='27',solint=86400.,gaincurve=True)
bandpass('ngc1333_1.ms','n1333.2.bcal',msselect='FIELD_ID==15 && DATA_DESC_ID==1',opacity=False,gaintable='n1333.2.gcal',refant='27',solint=86400.,gaincurve=True)
bptime = time.time()

###################################
### Transfer the flux density  from flux calibrater to gain calibraters  
print '--Fluxscale--'
default('fluxscale')
fluxscale('ngc1333_1.ms',caltable='n1333.1.gcal',fluxtable='n1333.1.fcal',reference='0542+498_1',transfer=['0336+323_1', '0319+415_1'])
fluxscale('ngc1333_1.ms',caltable='n1333.2.gcal',fluxtable='n1333.2.fcal',reference='0542+498_2',transfer=['0336+323_2', '0319+415_2'])
fstime = time.time()

###################################
###correct the visibilities of target with the previous solutions derived
###including the gain curve correction
print '--Correct--'
default('correct')
correct(vis='ngc1333_1.ms',msselect='FIELD_ID IN [0,1,2,3,4,5] && DATA_DESC_ID == 0',gaincurve=True,opacity=False,gaintable='n1333.1.fcal',gainselect='FIELD_ID IN [0]',bptable='n1333.1.bcal',calwt=False)
correct(vis='ngc1333_1.ms',msselect='FIELD_ID IN [6,7,8,9,10,11] && DATA_DESC_ID == 1',gaincurve=True,opacity=False,gaintable='n1333.2.fcal',gainselect='FIELD_ID IN [6]',bptable='n1333.2.bcal',calwt=False)
correcttime = time.time()


print '--Split (target and cals) --'
default('split')
######split out the corrected data of target 
split('ngc1333_1.ms','src.t2may.ms',fieldid=[1,2,3,4,5,7,8,9,10,11],spwid=-1,nchan=63,start=0,step=1,datacolumn='corrected')
##### split out calibraters
split('ngc1333_1.ms','gcal.t1.2may.ms',fieldid=0,spwid=0,nchan=63,start=0,step=1,datacolumn='corrected')
split('ngc1333_1.ms','gcal.t2.2may.ms',fieldid=6,spwid=1,nchan=63,start=0,step=1,datacolumn='corrected')
splitcaltime = time.time()

print '*** 08 MAY ***'
print '--Import--'
default('importvla')
###################################################
## Fill data from UVFITS format to measurementset
importuvfits(fitsfile=datapath+'N1333_2.UVFITS',vis='ngc1333_2.ms')
importtime2 = time.time()
#print '--Observation summary--'
#listobs('ngc1333_2.ms')
#listtime = time.time()
print '--Flag Autocorrelations--'
######################
###flag autocorrelations
default('flagautocorr')
flagautocorr('ngc1333_2.ms')
####flag all data from antennas 8,9,18
default('flagdata')
flagdata('ngc1333_2.ms',antennaid=[8,9,18],fieldid=-1)
####flag antenna 14 for timerange
default('flagdata')
flagdata('ngc1333_2.ms',antennaid=14,timerange=['08-MAY-2003/00:00:00','08-MAY-2003/20:00:00'],fieldid=-1)
######flag the last data channels
default('flagdata')
flagdata('ngc1333_2.ms',chans=[59,60,61,62],fieldid=-1,spwid=-1)
flagtime2=time.time()

print '--Setjy--'
default('setjy')
setjy('ngc1333_2.ms',field='0542+498_1')
setjy('ngc1333_2.ms',field='0542+498_2')
setjytime2 = time.time()

print '--Gaincal--'
default('gaincal')
gaincal('ngc1333_2.ms','n1333b.1.gcal',mode='channel',nchan=55,start=4,msselect='FIELD_ID IN [0,12,14] && DATA_DESC_ID IN [0]',opacity=False,solint=0.,refant='27')
gaincal('ngc1333_2.ms','n1333b.2.gcal',mode='channel',nchan=55,start=4,msselect='FIELD_ID IN [6,13,15] && DATA_DESC_ID IN [1]',opacity=False,solint=0.,refant='27')
gaintime2 = time.time()

print '--Bandpass--'
default('bandpass')
bandpass('ngc1333_2.ms','n1333b.1.bcal',msselect='FIELD_ID==14 && DATA_DESC_ID==0',opacity=False,gaintable='n1333b.1.gcal',refant='27',nchan=0,start=0,step=1,solint=86400.)
bandpass('ngc1333_2.ms','n1333b.2.bcal',msselect='FIELD_ID==15 && DATA_DESC_ID==1',opacity=False,gaintable='n1333b.2.gcal',refant='27',nchan=0,start=0,step=1,solint=86400.)
bptime2 = time.time()

print '--Fluxscale--'
default('fluxscale')
fluxscale('ngc1333_2.ms',caltable='n1333b.1.gcal',fluxtable='n1333b.1.fcal',reference='0542+498_1',transfer=['0336+323_1', '0319+415_1'])
fluxscale('ngc1333_2.ms',caltable='n1333b.2.gcal',fluxtable='n1333b.2.fcal',reference='0542+498_2',transfer=['0336+323_2', '0319+415_2'])
fstime2 = time.time()

#### Do the correction for the above solutions of bandpass and gain,
#### including gain curve
print '--Correct--'
default('correct')
correct(vis='ngc1333_2.ms',msselect='FIELD_ID IN [0,1,2,3,4,5] && DATA_DESC_ID == 0',gaincurve=True,opacity=False,gaintable='n1333b.1.fcal',gainselect='FIELD_ID IN [0]',bptable='n1333b.1.bcal',calwt=False)
correct(vis='ngc1333_2.ms',msselect='FIELD_ID IN [6,7,8,9,10,11] && DATA_DESC_ID == 1',gaincurve=True,opacity=False,gaintable='n1333b.2.fcal',gainselect='FIELD_ID IN [6]',bptable='n1333b.2.bcal',calwt=False)
correcttime2 = time.time()

###split corrected target and calibraters 
print '--Split (target and cals) --'
default('split')
split('ngc1333_2.ms','src.t8may.ms',fieldid=[1,2,3,4,5,7,8,9,10,11],spwid=-1,nchan=63,start=0,step=1,datacolumn='corrected')
split('ngc1333_2.ms','gcal.t1.8may.ms',fieldid=0,spwid=0,nchan=63,start=0,step=1,datacolumn='corrected')
split('ngc1333_2.ms','gcal.t2.8may.ms',fieldid=6,spwid=1,nchan=63,start=0,step=1,datacolumn='corrected')
splitcaltime2 = time.time()

#### merge the 2 corrected data sets
print '--Concatenate data sets--'
os.system('cp -r src.t2may.ms ./ngc1333_both.ms')
ms.open('ngc1333_both.ms',nomodify=False)
ms.concatenate('src.t8may.ms',freqtol='10MHz',dirtol='1arcsec')
ms.close()

#default('flagdata')
#flagdata('ngc1333_both.ms',chans=[59,60,61,62],fieldid=-1,spwid=-1)
#default('flagdata')
#flagdata('ngc1333_both.ms',timerange=['02-MAY-2003/21:50:40','02-MAY-2003/22:01:10'],fieldid=-1)
#default('flagdata')
#flagdata('ngc1333_both.ms',antennaid=13,timerange=['02-MAY-2003/18:50:50','02-MAY-2003/19:13:30'],fieldid=-1)
concattime=time.time()

#####Imaging using the imager tool
print '--Image data--'
im.open('ngc1333_both.ms')
###selecting data from 10 fields 0-9 all channels of the 2 spectral windows
im.selectvis(field=range(0,10),spwid=[0,1],nchan=63,start=0,step=1)
###image centered of field 0, making 18 channels averaging 5 data channels to
### 1 image channel
im.defineimage(nx=800,ny=800,cellx='0.5arcsec',stokes='I',mode='channel',nchan=18,start=2,step=5,phasecenter=0,spwid=[0,1])
### natural weighting
im.weight(type='natural')
### Use default primary beam
im.setvp(dovp=True,usedefaultvp=True,dosquint=False)
### Use the newer mosaic gridder where mosaicing is done in convolution
im.setoptions(padding=1.0,ftmachine='mosaic')
#### Define up to which point in the beam image will be shown 
im.setmfcontrol(scaletype='NONE',minpb=0.07,constpb=0.4)
#### clean with niter of 1 ...basically doing a just a dirty mosaic
im.clean(algorithm='mfclark',niter=1,gain=0.1,threshold=5.,model='src.tall.cln.model',mask='',image='src.task.image',residual='src.tall.cln.resid',npercycle=30)
im.close()
imagetime=time.time()

# Do moments
print '--Calculate moments--'
ia.open('src.task.image')
ia.moments(outfile='src.tmom0.red',moments=0,axis=3,mask='indexin(4,[3:9])',includepix=[0.003,100.0])
ia.moments(outfile='src.tmom0.blu',moments=0,axis=3,mask='indexin(4,[10:16])',includepix=[0.003,100.0])
ia.moments(outfile='src.tmom0.all',moments=0,axis=3,mask='indexin(4,[3:16])',includepix=[0.003,100.0])
ia.moments(outfile='src.tmom1.all',moments=1,axis=3,mask='indexin(4,[3:16])',includepix=[0.02,100.0])
ia.close()
momenttime=time.time()

endProc = time.clock()
endTime = time.time()

# Regression

ms.open('src.t2may.ms')
src_2may=max(ms.range(["amplitude"]).get('amplitude'))
ms.close
ms.open('src.t8may.ms')
src_8may=max(ms.range(["amplitude"]).get('amplitude'))
ms.close
ms.open('gcal.t1.8may.ms')
gcal1_8may=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('gcal.t2.8may.ms')
gcal2_8may=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('gcal.t1.2may.ms')
gcal1_2may=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('gcal.t2.2may.ms')
gcal2_2may=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ia.open('src.tmom0.all')
statistics=ia.statistics()
thistest_immax=statistics['max'][0]
thistest_imrms=statistics['rms'][0]
# 7499.6601 (n1333_both.ms)
src2may=1010.402
src8may=273.024
cal1_2may=298.137
cal2_2may=253.087
cal1_8may=59.37
cal2_8may=281.055
immax=1.355
imrms=0.191

diff_cal1_2may=abs((cal1_2may-gcal1_2may)/cal1_2may)
diff_cal2_2may=abs((cal2_2may-gcal2_2may)/cal2_2may)
diff_cal1_8may=abs((cal1_8may-gcal1_8may)/cal1_8may)
diff_cal2_8may=abs((cal2_8may-gcal2_8may)/cal2_8may)
diff_src_2may=abs((src2may-src_2may)/src2may)
diff_src_8may=abs((src8may-src_8may)/src8may)
diff_immax=abs((immax-thistest_immax)/immax)
diff_imrms=abs((imrms-thistest_imrms)/imrms)

import datetime
datestring=datetime.datetime.isoformat(datetime.datetime.today())
outfile='ngc1333.'+datestring+'.log'
logfile=open(outfile,'w')

print >>logfile,''
print >>logfile,'********** Data Summary *********'
print >>logfile,'*********************************'
print >>logfile,''
print >>logfile,'********** Regression ***********'
print >>logfile,'*                               *'
if (diff_cal1_2may < 0.05): print >>logfile,'* Passed cal1 max amplitude test (2may) *'
print >>logfile,'* Cal1 max amp (2may)'+str(gcal1_2may)
if (diff_cal2_2may < 0.05): print >>logfile,'* Passed cal2 max amplitude test (2may) *'
print >>logfile,'* Cal2 max amp (2may)'+str(gcal2_2may)
if (diff_cal1_8may < 0.05): print >>logfile,'* Passed cal1 max amplitude test (8may) *'
print >>logfile,'* Cal1 max amp (8may) '+str(gcal1_8may)
if (diff_cal2_8may < 0.05): print >>logfile,'* Passed cal2 max amplitude test (8may) *'
print >>logfile,'* Cal2 max amp (8may) '+str(gcal2_8may)
if (diff_src_2may < 0.05): print >>logfile,'* Passed src max amplitude test (2may) *'
print >>logfile,'* Src max amp (2may) '+str(src_2may)
if (diff_src_8may < 0.05): print >>logfile,'* Passed src max amplitude test (8may) *'
print >>logfile,'* Src max amp (8may) '+str(src_8may)
print >>logfile,'* Image max '+str(thistest_immax)
if (diff_imrms < 0.05): print >>logfile,'* Passed image rms test         *'
print >>logfile,'* Image rms '+str(thistest_imrms)
if ((diff_cal1_2may<0.05) & (diff_cal2_2may<0.05) & (diff_cal1_8may<0.05) & (diff_cal2_8may<0.05) & (diff_src_2may<0.05) & (diff_src_8may<0.05) & (diff_immax<0.05) & (diff_imrms<0.05)): 
	regstate=True
	print >>logfile,'---'
	print >>logfile,'Passed Regression test for NGC1333'
	print >>logfile,'---'
else: 
	regstate=False
	print >>logfile,'----FAILED Regression test for NGC1333'
print >>logfile,'*********************************'

print >>logfile,''
print >>logfile,''
print >>logfile,'********* Benchmarking *****************'
print >>logfile,'*                                      *'
print >>logfile,'Total wall clock time was: '+str(endTime - startTime)
print >>logfile,'Total CPU        time was: '+str(endProc - startProc)
print >>logfile,'Processing rate MB/s  was: '+str(240.3/(endTime - startTime))
print >>logfile,'* Breakdown:                           *'
print >>logfile,'*   import       time was: '+str(importtime-startTime)
print >>logfile,'*   flagdata     time was: '+str(flagtime-importtime)
print >>logfile,'*   setjy        time was: '+str(setjytime-flagtime)
print >>logfile,'*   gaincal      time was: '+str(gaintime-setjytime)
print >>logfile,'*   bandpass     time was: '+str(bptime-gaintime)
print >>logfile,'*   fluxscale    time was: '+str(fstime-bptime)
print >>logfile,'*   correct      time was: '+str(correcttime-fstime)
print >>logfile,'*   split        time was: '+str(splitcaltime-correcttime)
print >>logfile,'*   import       time was: '+str(importtime2-splitcaltime)
print >>logfile,'*   flagdata     time was: '+str(flagtime2-importtime2)
print >>logfile,'*   setjy        time was: '+str(setjytime2-flagtime2)
print >>logfile,'*   gaincal      time was: '+str(gaintime2-setjytime2)
print >>logfile,'*   bandpass     time was: '+str(bptime2-gaintime2)
print >>logfile,'*   fluxscale    time was: '+str(fstime2-bptime2)
print >>logfile,'*   correct      time was: '+str(correcttime2-fstime2)
print >>logfile,'*   split        time was: '+str(splitcaltime2-correcttime2)
print >>logfile,'*   concatenate  time was: '+str(concattime-splitcaltime2)
print >>logfile,'*   image        time was: '+str(imagetime-concattime)
print >>logfile,'*   moments      time was: '+str(momenttime-imagetime)
print >>logfile,'*****************************************'
#
logfile.close()