|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################### 00002 # # 00003 # Regression/Benchmarking Script for NGC 7538 # 00004 # # 00005 ############################################### 00006 00007 import time 00008 import os 00009 00010 os.system('rm -rf ngc7538d* 1328.* 2229.cont2* ap314.* ngc7538*.ms') 00011 00012 startTime = time.time() 00013 startProc = time.clock() 00014 00015 datapath=os.environ.get('CASAPATH').split()[0]+'/data/regression/ATST1/NGC7538/' 00016 00017 print '--Import--' 00018 default('importvla') 00019 importvla(archivefiles=[datapath+'AP314_A950519.xp1', 00020 datapath+'AP314_A950519.xp2', 00021 datapath+'AP314_A950519.xp3'], 00022 vis='ngc7538.ms', bandname='K', frequencytol=10000000.0) 00023 importtime = time.time() 00024 print '--Observation summary--' 00025 listobs(vis='ngc7538.ms') 00026 #listtime = time.time() 00027 print '--Flag auto-correlations--' 00028 default('tflagdata') 00029 tflagdata(vis='ngc7538.ms', mode='manual', autocorr=True) 00030 flagtime = time.time() 00031 print '--Setjy--' 00032 default('setjy') 00033 setjy(vis='ngc7538.ms',field='0',standard='Perley-Taylor 99',scalebychan=False) #set flux density for 1331+305 (3C286) 00034 setjytime = time.time() 00035 print '--Gaincal--' 00036 default('gaincal') 00037 gaincal(vis='ngc7538.ms', caltable='ap314.gcal', 00038 field='<2', spw='0~1:2~56', gaintype='G', 00039 opacity=0.08,solint='inf', combine='', refant='VA19') 00040 gaintime = time.time() 00041 print '--Bandpass--' 00042 default('bandpass') 00043 bandpass(vis='ngc7538.ms', caltable='1328.bcal', 00044 field='0', opacity=0.08, 00045 gaintable='ap314.gcal', interp='nearest', 00046 refant='VA19') 00047 bptime = time.time() 00048 print '--Fluxscale--' 00049 default('fluxscale') 00050 fluxscale(vis='ngc7538.ms', caltable='ap314.gcal', fluxtable='ap314.fluxcal', 00051 reference=['1328+307'], transfer=['2229+695']) 00052 fstime = time.time() 00053 print '--Apply Cal--' 00054 default('applycal') 00055 applycal(vis='ngc7538.ms', 00056 field='1~5', 00057 opacity=0.08, 00058 gaintable=['ap314.fluxcal', '1328.bcal'], 00059 gainfield='1') 00060 00061 correcttime = time.time() 00062 00063 print '--Split (fluxcal data)--' 00064 default('split') 00065 split(vis='ngc7538.ms', outputvis='ngc7538_cal.split.ms', 00066 field='0',spw='0:0~61', datacolumn='model') 00067 00068 print '--Split (continuum)--' 00069 default('split') 00070 # This _averages_ 00071 split(vis='ngc7538.ms', outputvis='ngc7538d.cont.ms', 00072 field='3',spw='0:2~56',width=[55], datacolumn='corrected') 00073 00074 print '--Split (mf cont,)--' 00075 default('split') 00076 split(vis='ngc7538.ms', outputvis='ngc7538.cont.ms', 00077 field='3,4,5',spw='0:2~56',width=[55], datacolumn='corrected') 00078 print '--Split (bandcal data)--' 00079 default('split') 00080 split(vis='ngc7538.ms', outputvis='2229.cont2.ms', 00081 field='1', spw='0:2~56,1:2~56',width=[55,55], datacolumn='corrected') 00082 splitcaltime = time.time() 00083 default('split') 00084 split(vis='ngc7538.ms',outputvis='ngc7538d.line.ms', 00085 field='3',spw='0:2~56,1:2~56',datacolumn='corrected') 00086 splitsrctime = time.time() 00087 00088 print '--Clean cal--' 00089 default('clean') 00090 clean(vis='2229.cont2.ms',imagename='2229.cont2',mode='channel', 00091 psfmode='hogbom',niter=6000,gain=0.1,threshold=8.,mask='', 00092 nchan=1,start=0,width=1,field='0',spw='0', 00093 imsize=[256,256],cell=[0.5,0.5], 00094 weighting='briggs',robust=0.5,imagermode='') 00095 cleantime1 = time.time() 00096 print '--Clean src cont--' 00097 default('clean') 00098 clean(vis='ngc7538d.cont.ms',imagename='ngc7538d.cont',mode='channel', 00099 psfmode='hogbom',niter=5000,gain=0.1,threshold=3.,mask='', 00100 nchan=1,start=0,width=1,field='0',spw='0', 00101 imsize=[1024,1024],cell=[0.5,0.5], 00102 weighting='briggs',robust=0.5,imagermode='') 00103 cleantime2 = time.time() 00104 print '--Clean src line--' 00105 default('clean') 00106 clean(vis='ngc7538d.line.ms',imagename='ngc7538d.cube',mode='channel', 00107 psfmode='hogbom',niter=5000,gain=0.1,threshold=30.,mask='', 00108 nchan=48,start=2,width=1,field='0',spw='0', 00109 imsize=[128,128],cell=[4.,4.], 00110 weighting='briggs',robust=2.,imagermode='', 00111 uvtaper=True, outertaper=['12.0arcsec','12.0arcsec', '0deg']) 00112 cleantime3 = time.time() 00113 # -- Not done in old regression but should be 00114 #print '--Contsub (image plane)--' 00115 #ia.open('ngc7538d.cube.image') 00116 #myim=ia.continuumsub('ngc7538d_subed.line.im','ngc7538d_res.cont.im',channels=range(0,48),fitorder=1) 00117 #ia.close() 00118 #myim.close() 00119 #contsubtime=time.time() 00120 #print '--View image--' 00121 #viewer('ngc5921_task.image') 00122 00123 endProc = time.clock() 00124 endTime = time.time() 00125 00126 # Regression 00127 00128 test_name_cal = 'NGC7538--Calibrater maximum amplitude test' 00129 test_name_src = 'NGC7538--Source maximum amplitude test' 00130 test_name_immax = 'NGC7538--Image maximum test' 00131 test_name_imrms = 'NGC7538--Image rms test' 00132 00133 ms.open('ngc7538_cal.split.ms') 00134 thistest_cal=max(ms.range(["amplitude"]).get('amplitude')) 00135 ms.close() 00136 ms.open('ngc7538d.line.ms') 00137 thistest_src=max(ms.range(["amplitude"]).get('amplitude')) 00138 ms.close() 00139 ia.open('ngc7538d.cube.image') 00140 # ia.statistics returns dictionary with 'return','statsout' 00141 # get the second value in the dictionary (statsout) 00142 statistics=ia.statistics() 00143 # note thistest_immax will be a list with one value 00144 thistest_immax=statistics['max'][0] 00145 # note thistest_imrms will be a list with one value 00146 thistest_imrms=statistics['rms'][0] 00147 00148 cal_max=2.413 00149 src_max=18.3638 00150 im_max=0.2606 00151 im_rms=0.0124 00152 00153 00154 diff_cal=abs((cal_max-thistest_cal)/cal_max) 00155 diff_src=abs((src_max-thistest_src)/src_max) 00156 diff_immax=abs((im_max-thistest_immax)/im_max) 00157 diff_imrms=abs((im_rms-thistest_imrms)/im_rms) 00158 00159 import datetime 00160 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00161 outfile='ngc7538.'+datestring+'.log' 00162 logfile=open(outfile,'w') 00163 00164 print >>logfile,'' 00165 print >>logfile,'********** Data Summary *********' 00166 print >>logfile,'* Observer: unavailable Project: AP314 *' 00167 print >>logfile,'* Observation: VLA(27 antennas) *' 00168 print >>logfile,'* Telescope Observation Date Observer Project *' 00169 print >>logfile,'* VLA [ 4.30759e+09, 4.30759e+09]unavailable AP314 *' 00170 print >>logfile,'* VLA [ 4.30759e+09, 4.30762e+09]unavailable AP314 *' 00171 print >>logfile,'* VLA [ 4.30762e+09, 4.30763e+09]unavailable AP314 *' 00172 print >>logfile,'*Data records: 838404 Total integration time = 36000 seconds *' 00173 print >>logfile,'* Observed from 09:23:45 to 19:23:45 *' 00174 print >>logfile,'*Fields: 6 *' 00175 print >>logfile,'* ID Name Right Ascension Declination Epoch *' 00176 print >>logfile,'* 0 1328+307 13:31:08.29 +30.30.33.04 J2000 *' 00177 print >>logfile,'* 1 2229+695 22:30:36.48 +69.46.28.00 J2000 *' 00178 print >>logfile,'* 2 NGC7538C 23:14:02.48 +61.27.14.86 J2000 *' 00179 print >>logfile,'* 3 NGC7538D 23:13:43.82 +61.27.00.18 J2000 *' 00180 print >>logfile,'* 4 NGC7538E 23:13:34.64 +61.27.26.44 J2000 *' 00181 print >>logfile,'* 5 NGC7538F 23:13:35.76 +61.28.33.66 J2000 *' 00182 print >>logfile,'* Data descriptions: 2 (2 spectral windows and 2 polarization setups) *' 00183 print >>logfile,'* ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs *' 00184 print >>logfile,'* 0 63 TOPO 23691.4682 118.164062 6152.34375 23694.4955 RR *' 00185 print >>logfile,'* 1 63 TOPO 23719.6063 118.164062 6152.34375 23722.6336 LL *' 00186 print >>logfile,'*********************************' 00187 print >>logfile,'' 00188 print >>logfile,'********** Regression ***********' 00189 print >>logfile,'* *' 00190 passfail = {True: '* Passed', 00191 False: '* FAILED'} 00192 print >> logfile, passfail[diff_cal < 0.05], 'cal max amplitude test *' 00193 print >>logfile,'* Cal max amp '+str(thistest_cal) 00194 print >> logfile, passfail[diff_src < 0.05], 'src max amplitude test *' 00195 print >>logfile,'* Src max amp '+str(thistest_src) 00196 print >> logfile, passfail[diff_immax < 0.05], 'image max test *' 00197 print >>logfile,'* Image max '+str(thistest_immax) 00198 print >> logfile, passfail[diff_imrms < 0.05], 'image rms test *' 00199 print >>logfile,'* Image rms '+str(thistest_imrms) 00200 if ((diff_src<0.05) & (diff_cal<0.05) & (diff_immax<0.05) & (diff_imrms<0.05)): 00201 regstate=True 00202 print >>logfile,'---' 00203 print >>logfile,'Passed Regression test for NGC7538' 00204 print >>logfile,'---' 00205 else: 00206 regstate=False 00207 print >>logfile,'----FAILED Regression test for NGC7538' 00208 print >>logfile,'*********************************' 00209 00210 print >>logfile,'' 00211 print >>logfile,'' 00212 print >>logfile,'********* Benchmarking *****************' 00213 print >>logfile,'* *' 00214 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00215 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00216 print >>logfile,'Processing rate MB/s was: '+str(240.3/(endTime - startTime)) 00217 print >>logfile,'* Breakdown: *' 00218 print >>logfile,'* import time was: '+str(importtime-startTime) 00219 print >>logfile,'* flagautocorr time was: '+str(flagtime-importtime) 00220 print >>logfile,'* setjy time was: '+str(setjytime-flagtime) 00221 print >>logfile,'* gaincal time was: '+str(gaintime-setjytime) 00222 print >>logfile,'* bandpass time was: '+str(bptime-gaintime) 00223 print >>logfile,'* fluxscale time was: '+str(fstime-bptime) 00224 print >>logfile,'* applycal time was: '+str(correcttime-fstime) 00225 print >>logfile,'* split-cal time was: '+str(splitcaltime-correcttime) 00226 print >>logfile,'* split-src time was: '+str(splitsrctime-splitcaltime) 00227 print >>logfile,'* clean-cal time was: '+str(cleantime1-splitsrctime) 00228 print >>logfile,'* clean-src-c time was: '+str(cleantime2-cleantime1) 00229 print >>logfile,'* clean-src-l time was: '+str(cleantime3-cleantime2) 00230 #print '* contsub time was: ',contsubtime-cleantime3 00231 print >>logfile,'*****************************************' 00232 print >>logfile,'basho (test cpu) time was: 500 seconds' 00233 00234 logfile.close() 00235
1.8.0
