###############################################
#                                             #
# Regression/Benchmarking Script for NGC 5921 #
#                                             #
###############################################

import time
import os

os.system('rm -rf ngc5921.*cal* ngc5921_task.* ngc5921_*.split* ngc5921.ms cleanbox.mask')

#pathname=`/bin/echo $AIPSPATH | awk '{ print $1 }'`
pathname=os.environ.get('AIPSPATH')
datapath=pathname+'/data/demo/NGC5921.fits'

restore()

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

print '--Import--'
default('importuvfits')
importuvfits(datapath,'ngc5921.ms')
importtime = time.time() 
print '--Flagautocorr--'
default('flagautocorr')
flagautocorr('ngc5921.ms')
flagtime = time.time()
print '--Setjy--'
default('setjy')
setjy('ngc5921.ms',fieldid=0)
setjytime = time.time()
print '--Gaincal--'
default('gaincal')
gaincal('ngc5921.ms','ngc5921.gcal',mode='channel',nchan=55,start=2,msselect='FIELD_ID <= 1',tau=0.0001,solint=0.,refant='15')
gaintime = time.time()
print '--Bandpass--'
default('bandpass')
bandpass('ngc5921.ms','ngc5921.bcal',msselect='FIELD_ID==0',tau=0.0001,gaintable='ngc5921.gcal',refant='15',nchan=0,start=0,step=1,solint=86400.0)
bptime = time.time()
print '--Fluxscale--'
default('fluxscale')
fluxscale('ngc5921.ms',caltable='ngc5921.gcal',fluxtable='ngc5921.fluxscale',reference=['1331+30500002_0'],transfer=['1445+09900002_0'])
fstime = time.time()
print '--Correct--'
default('correct')
correct('ngc5921.ms',msselect='FIELD_ID IN [1,2]',tau=0.0001,gaintable='ngc5921.fluxscale',gainselect='FIELD_ID==1',bptable='ngc5921.bcal')
correcttime = time.time()
print '--Split (cal data)--'
default('split')
split('ngc5921.ms','ngc5921_cal.split.ms',fieldid=1,spwid=0,nchan=63,start=0,step=1,datacolumn='CORRECTED_DATA')
splitcaltime = time.time()
print '--Continuum Subtract--'
default('contsub')
contsub('ngc5921.ms',fieldid=2,spwid=0,channels=range(4,7)+range(50,60),solint=0.0,fitmode='subtract')
contsubtime = time.time()
print '--Split (src data)--'
default('split')
split('ngc5921.ms','ngc5921_src.split.ms',fieldid=2,spwid=0,nchan=46,start=5,step=1,datacolumn='CORRECTED_DATA')
splitsrctime = time.time()
print '--Clean--'
default('clean')
clean('ngc5921_src.split.ms','ngc5921_task','channel','hogbom',niter=6000,gain=0.1,threshold=8.,mask='',nchan=46,start=0,step=1,fieldid=0,imsize=[256,256],cell=[15.,15.],weighting='briggs',rmode='norm',robust=0.5)
#clean('ngc5921_src.split.ms','ngc5921_task','channel','hogbom',niter=6000,gain=0.1,threshold=8.,mask='',nchan=46,start=0,step=1,fieldid=0,imsize=[256,256],cell=[15.,15.],weighting='briggs',rmode='norm',robust=0.5,cleanbox='regionfile.txt')
cleantime = time.time()


#print '--View image--'
#viewer('ngc5921_task.image')

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

# Regression

test_name_cal = 'NGC5921--Calibrater maximum amplitude test'
test_name_src = 'NGC5921--Source maximum amplitude test'
test_name_immax = 'NGC5921--Image maximum test'
test_name_imrms = 'NGC5921--Image rms test'

ms.open('ngc5921_cal.split.ms')
thistest_cal=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('ngc5921_src.split.ms')
thistest_src=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ia.open('ngc5921_task.image')
statistics=ia.statistics()
thistest_immax=statistics['max'][0]
# note thistest_imrms will be a list with one value 
thistest_imrms=statistics['rms'][0]

#cal max amplitude = 34.0338668823
#src max amplitude =  1.37963354588
#im  max           = 0.052414759993553162
#im  rms           = 0.0020218724384903908
#data set size     = 35.1 MB UVFITS file

diff_cal=abs((34.0338668823-thistest_cal)/34.0338668823)
diff_src=abs((1.37963354588-thistest_src)/1.37963354588)
diff_immax=abs((0.052414759993553162-thistest_immax)/0.052414759993553162)
diff_imrms=abs((0.0020218724384903908-thistest_imrms)/0.0020218724384903908)

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

print >>logfile,''
print >>logfile,''
print >>logfile,'********** Data Summary *********'
print >>logfile,'*   Observer: TEST     Project:                                             *'
print >>logfile,'* Observation: VLA(28 antennas)                                             *'
print >>logfile,'* Data records: 22653  Total integration time = 5280 seconds                *'
print >>logfile,'* Observed from   09:19:00   to   10:47:00                                  *'
print >>logfile,'* Fields: 3                                                                 *'
print >>logfile,'*  ID   Name          Right Ascension  Declination   Epoch                  *'
print >>logfile,'*  0    1331+30500002_013:31:08.29      +30.30.32.96  J2000                 *'
print >>logfile,'*  1    1445+09900002_014:45:16.47      +09.58.36.07  J2000                 *'
print >>logfile,'*  2    N5921_2       15:22:00.00      +05.04.00.00  J2000                  *'
print >>logfile,'* Data descriptions: 1 (1 spectral windows and 1 polarization setups)       *'
print >>logfile,'*  ID  #Chans Frame Ch1(MHz)    Resoln(kHz) TotBW(kHz)  Ref(MHz)    Corrs   *'  
print >>logfile,'*  0       63 LSRK  1412.68608  24.4140625  1550.19688  1413.44902  RR  LL  *'
print >>logfile,'*********************************'
print >>logfile,''
print >>logfile,'********** Regression ***********'
print >>logfile,'*                               *'
if (diff_cal < 0.05): print >>logfile,'* Passed cal max amplitude test '
print >>logfile,'*  Cal max amp '+str(thistest_cal)
if (diff_src < 0.05): print >>logfile,'* Passed src max amplitude test '
print >>logfile,'*  Src max amp '+str(thistest_src)
if (diff_immax < 0.05): print >>logfile,'* Passed image max test'
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_src<0.05) & (diff_cal<0.05) & (diff_immax<0.05) & (diff_imrms<0.05)): 
	regstate=True
	print >>logfile,'---'
	print >>logfile,'Passed Regression test for NGC5921'
	print >>logfile,'---'
else: 
	regstate=False
	print >>logfile,'----FAILED Regression test for NGC5921'
print >>logfile,'*********************************'

print >>logfile,''
print >>logfile,''
print >>logfile,'********* Benchmarking *****************'
print >>logfile,'*                                      *'
print >>logfile,'Total wall clock time was: ', endTime - startTime
print >>logfile,'Total CPU        time was: ', endProc - startProc
print >>logfile,'Processing rate MB/s  was: ', 35.1/(endTime - startTime)
print >>logfile,'* Breakdown:                           *'
print >>logfile,'*   import       time was: '+str(importtime-startTime)
print >>logfile,'*   flagautocorr 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-cal    time was: '+str(splitcaltime-correcttime)
print >>logfile,'*   contsub      time was: '+str(contsubtime-splitcaltime)
print >>logfile,'*   split-src    time was: '+str(splitsrctime-contsubtime)
print >>logfile,'*   clean        time was: '+str(cleantime-splitsrctime)
print >>logfile,'*****************************************'
print >>logfile,'basho (test cpu) time was: 55 seconds'

logfile.close()