|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################### 00002 # # 00003 # Regression/Benchmarking Script for B0319 # 00004 # # 00005 ############################################### 00006 00007 import time 00008 import os 00009 00010 os.system('rm -rf B0319_0317.ms B0319.K B0319.Mt B0319.MFt B0319*.png n1333.ms') 00011 00012 pathname=os.environ.get('CASAPATH').split()[0] 00013 datapath = os.environ.get('CASAPATH').split()[0]+'/data/regression/ATST4/B0319/N1333_1.UVFITS' 00014 00015 startTime = time.time() 00016 startProc = time.clock() 00017 00018 00019 # Baseline-based calibration of N1333 calibrater 00020 # VLA baseline 3-17 02-May-2003 00021 # MeasurementSet Name: B0319_0317.ms MS Version 2 00022 # 00023 # Observer: AW602 Project: 00024 #Observation: VLA 00025 #Data records: 62 Total integration time = 609.999 seconds 00026 # Observed from 02-May-2003/20:09:40 to 02-May-2003/20:19:50 00027 # 00028 # ObservationID = 1 ArrayID = 1 00029 # Date Timerange Scan FldId FieldName DataDescIds 00030 # 02-May-2003/20:09:40.0 - 20:19:50.0 49 1 0319+415_1 [1] 00031 #Fields: 1 00032 # ID Name Right Ascension Declination Epoch 00033 # 1 0319+415_1 03:19:48.16 +41.30.42.10 J2000 00034 #Data descriptions: 1 (1 spectral windows and 1 polarization setups) 00035 # ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 00036 # 1 63 LSRK 43416.2392 97.65625 6250 43419.2666 RR LL 00037 #Feeds: 29: printing first row only 00038 # Antenna Spectral Window # Receptors Polarizations 00039 # 1 -1 2 [ R, L] 00040 #Antennas: 2: 00041 # ID Name Station Diam. Long. Lat. 00042 # 3 3 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1 00043 # 17 17 VLA:E3 25.0 m -107.37.02.8 +33.54.00.5 00044 00045 # Examine observed data 00046 00047 00048 # N1333 00049 # ===== 00050 #FIELDID NAME PURPOSE COMMENT 00051 #------- ---- ------- ------- 00052 #1 0336_323_1 Gain Calibrater 00053 #13 0542+498_1 Flux Calibrater (3C147) 00054 #15 0319+415_1 Band Pass Calibrater (3C84) 00055 00056 startTime = time.time() 00057 startProc = time.clock() 00058 00059 print '--Import--' 00060 default('importuvfits') 00061 importuvfits(fitsfile=datapath,vis='n1333.ms',antnamescheme='new') 00062 importtime=time.time() 00063 print '--Split Data--' 00064 default('split') 00065 split(vis='n1333.ms',outputvis='B0319_0317.ms',datacolumn='data',field='14',antenna='VA03 & VA17') 00066 splittime=time.time() 00067 00068 print '--Plot antenna array and uv coverage--' 00069 default('plotxy') 00070 plotxy(vis='B0319_0317.ms',xaxis='x',subplot=121,markersize=8) 00071 plotxy(vis='B0319_0317.ms',xaxis='u',subplot=122,markersize=5) 00072 plotanttime=time.time() 00073 00074 print '--Plot visibility phase/amp versus uv distance--' 00075 default('plotxy') 00076 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp', 00077 datacolumn='data',subplot=211,clearpanel='All') 00078 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='phase', 00079 datacolumn='data',subplot=212,clearpanel='Auto') 00080 plotuvdisttime=time.time() 00081 00082 default('plotxy') 00083 plotxy(vis='B0319_0317.ms',xaxis='time',yaxis='amp', 00084 datacolumn='data',subplot=211,clearpanel='All') 00085 plotxy(vis='B0319_0317.ms',xaxis='time',yaxis='phase', 00086 datacolumn='data',subplot=212,clearpanel='Auto') 00087 plottimetime=time.time() 00088 00089 default('plotxy') 00090 plotxy(vis='B0319_0317.ms',xaxis='channel',yaxis='amp', 00091 datacolumn='data',subplot=211,clearpanel='All') 00092 plotxy(vis='B0319_0317.ms',xaxis='channel',yaxis='phase', 00093 datacolumn='data',subplot=212,clearpanel='Auto') 00094 plotchanneltime=time.time() 00095 00096 #Example data flagging 00097 #Non-interactive flagging has been removed from the plotxy task. 00098 #So this test as been temporarily removed. 00099 #default('flagxy') 00100 #plotxy(vis='B0319_0317.ms',xaxis='channel',yaxis='amp',datacolumn='data', 00101 # flagregion=[58,65,0,2.],enableflag=True,flagmode='f',clearpanel=True) 00102 #plotxy(vis='B0319_0317.ms',xaxis='channel',yaxis='amp',datacolumn='data', 00103 # flagregion=[0,3,0,2.],enableflag=True,flagmode='f',clearpanel=False) 00104 plotflagtime=time.time() 00105 00106 #Example spectral averaging 00107 default('plotxy') 00108 plotxy(vis='B0319_0317.ms',xaxis='time',yaxis='amp',datacolumn='data', 00109 subplot=111,clearpanel='All') 00110 plotxy(vis='B0319_0317.ms',xaxis='time',yaxis='amp',datacolumn='data', 00111 subplot=111,spw='0~1:0~60^20',overplot=True, 00112 plotsymbol='bo',clearpanel='All') 00113 plotavertime=time.time() 00114 00115 #Calibrate data 00116 clearcal(vis='B0319_0317.ms') 00117 default('blcal') 00118 blcal(vis='B0319_0317.ms',caltable='B0319.Mt', 00119 solint='3s',combine='',gaincurve=False,opacity=0.0) 00120 default('blcal') 00121 blcal(vis='B0319_0317.ms',caltable='B0319.MFt', 00122 gaintable='B0319.Mt',interp='nearest', 00123 solint='inf',combine='scan', 00124 gaincurve=False,opacity=0.0,freqdep=True) 00125 default('applycal') 00126 applycal(vis='B0319_0317.ms', 00127 gaintable=['B0319.Mt','B0319.MFt'], 00128 gaincurve=False,opacity=0.0) 00129 calibratetime=time.time() 00130 00131 #Compare observed 'data' and 'corrected' data 00132 default('plotxy') 00133 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp',datacolumn='data', 00134 plotsymbol=',',plotcolor='blue',subplot=121,clearpanel='All') 00135 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp',datacolumn='corrected', 00136 subplot=121,plotsymbol='+',plotcolor='red',overplot=True,clearpanel='Auto') 00137 #spectral average comparison, channel selection done with SPW selection 00138 default('plotxy') 00139 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp',datacolumn='data', 00140 subplot=122,plotsymbol=',',plotcolor='blue',spw='0~1:5~54^50') 00141 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp',datacolumn='corrected', 00142 subplot=122,plotsymbol='+',plotcolor='red',spw='0~1:5~54^50', 00143 overplot=True) 00144 plotcomparetime=time.time() 00145 00146 ##Compare M,MF with K 00147 #plotxy('B0319_0317.ms','time','phase','data',nchan=1,start=5,width=50,plotsymbol='ro') 00148 #plotxy('B0319_0317.ms','time','phase','corrected',nchan=1,start=5,width=50,plotsymbol='bo',overplot=True) 00149 ##Derive K 00150 #fringecal('B0319_0317.ms','B0319.K',solint=0.,gaincurve=False,opacity=0.0) 00151 #correct('B0319_0317.ms',gaintable='B0319.K') 00152 #plotxy('B0319_0317.ms','time','phase','corrected',nchan=1,start=5,width=50,plotsymbol='go',overplot=True) 00153 00154 # M,MF is better in this case - re-correct the data 00155 #clearcal('B0319_0317.ms') 00156 #blcal('B0319_0317.ms',caltable='B0319.Mt',solint=3.,gaincurve=False,opacity=0.0) 00157 #blcal('B0319_0317.ms',caltable='B0319.MFt',solint=30000.,gaincurve=False,opacity=0.0,freqdep=True) 00158 #correct('B0319_0317.ms',gaintable=['B0319.Mt','B0319.MFt']) 00159 00160 #Examine the calibration solutions 00161 default('plotcal') 00162 plotcal(caltable='B0319.Mt',yaxis='phase',subplot=121,plotsymbol='bo',clearpanel='All') 00163 plotcal(caltable='B0319.Mt',yaxis='amp',subplot=122,plotsymbol='bo', 00164 overplot=True,clearpanel='Auto') 00165 00166 default('plotcal') 00167 plotcal(caltable='B0319.MFt',yaxis='phase', 00168 subplot=121,plotsymbol='bo',clearpanel='All') 00169 plotcal(caltable='B0319.MFt',yaxis='amp', 00170 subplot=122,plotsymbol='bo',overplot=True,clearpanel='Auto') 00171 plotcaltime=time.time() 00172 00173 #default('plotcal') 00174 #plotcal('B0319.K','delay',subplot=121,plotsymbol='go',clearpanel=True) 00175 #plotcal('B0319.K','delayrate',subplot=122,plotsymbol='go',clearpanel=False) 00176 00177 #Examine corrected data 00178 default('plotxy') 00179 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp',datacolumn='corrected', 00180 subplot=121,clearpanel='All') 00181 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='phase',datacolumn='corrected', 00182 subplot=122,clearpanel='Auto') 00183 00184 default('plotxy') 00185 plotxy(vis='B0319_0317.ms',xaxis='time',yaxis='amp',datacolumn='corrected', 00186 subplot=121,clearpanel='All') 00187 plotxy(vis='B0319_0317.ms',xaxis='time',yaxis='phase',datacolumn='corrected', 00188 subplot=122,clearpanel='Auto') 00189 00190 default('plotxy') 00191 plotxy(vis='B0319_0317.ms',xaxis='channel',yaxis='amp',datacolumn='corrected', 00192 subplot=121,clearpanel='All') 00193 plotxy(vis='B0319_0317.ms',xaxis='channel',yaxis='phase',datacolumn='corrected', 00194 subplot=122,clearpanel='Auto') 00195 plotcorrectedtime=time.time() 00196 00197 # uv model fit the data 00198 default('uvmodelfit') 00199 uvmodelfit(vis='B0319_0317.ms',niter=5,comptype='P', 00200 sourcepar=[0.5,.1,.1],outfile='test.cl') 00201 uvmodelfittime=time.time() 00202 00203 # now use component list to generate model data 00204 default('ft') 00205 ft(vis='B0319_0317.ms',complist='test.cl') 00206 fttime=time.time() 00207 00208 # Plot 00209 default('plotxy') 00210 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp',datacolumn='corrected', 00211 spw='0~1:5~54^50',plotsymbol='b,',clearpanel='All') 00212 plotxy(vis='B0319_0317.ms',xaxis='uvdist',yaxis='amp',datacolumn='model', 00213 spw='0~1:5~54^50',plotsymbol='ro',overplot=True,clearpanel='Auto') 00214 plotmodeltime=time.time() 00215 00216 endProc = time.clock() 00217 endTime = time.time() 00218 00219 #Regression 00220 00221 #test_name_amp = 'B0319 -- visibility max amplitude test' 00222 #test_name_ph = 'B0319 -- visibility max phase test' 00223 #test_name_mod = 'B0319 -- visibility max model test' 00224 00225 ms.open('B0319_0317.ms') 00226 #thistest_amp=pl.mean(ms.range(['corrected_amplitude']).get('corrected_amplitude')) 00227 #thistest_ph =pl.mean(ms.range(['corrected_phase']).get('corrected_phase')) 00228 thistest_mod=pl.mean(ms.range(['model_amplitude']).get('model_amplitude')) 00229 00230 #model amplitude 00231 model_amp=1.0 00232 00233 diff_mod=abs((model_amp-thistest_mod)/model_amp) 00234 00235 print '***' 00236 print 'model_amp ',model_amp 00237 print 'thistest_mod ',thistest_mod 00238 print '***' 00239 00240 import datetime 00241 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00242 outfile='B0319.'+datestring+'.log' 00243 logfile=open(outfile,'w') 00244 00245 print >>logfile,' Baseline-based calibration of N1333 calibrater' 00246 print >>logfile,' VLA baseline 3-17 02-May-2003' 00247 print >>logfile,' MeasurementSet Name: B0319_0317.ms MS Version 2' 00248 print >>logfile,'' 00249 print >>logfile,' Observer: AW602 Project:' 00250 print >>logfile,'Observation: VLA' 00251 print >>logfile,'Data records: 62 Total integration time = 609.999 seconds' 00252 print >>logfile,' Observed from 02-May-2003/20:09:40 to 02-May-2003/20:19:50' 00253 print >>logfile,'' 00254 print >>logfile,' ObservationID = 1 ArrayID = 1' 00255 print >>logfile,' Date Timerange Scan FldId FieldName DataDescIds' 00256 print >>logfile,' 02-May-2003/20:09:40.0 - 20:19:50.0 49 1 0319+415_1 [1]' 00257 print >>logfile,'Fields: 1' 00258 print >>logfile,' ID Name Right Ascension Declination Epoch' 00259 print >>logfile,' 1 0319+415_1 03:19:48.16 +41.30.42.10 J2000' 00260 print >>logfile,'Data descriptions: 1 (1 spectral windows and 1 polarization setups)' 00261 print >>logfile,' ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs' 00262 print >>logfile,' 1 63 LSRK 43416.2392 97.65625 6250 43419.2666 RR LL' 00263 print >>logfile,'Feeds: 29: printing first row only' 00264 print >>logfile,' Antenna Spectral Window # Receptors Polarizations' 00265 print >>logfile,' 1 -1 2 [ R, L]' 00266 print >>logfile,'Antennas: 2:' 00267 print >>logfile,' ID Name Station Diam. Long. Lat.' 00268 print >>logfile,' 3 3 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1' 00269 print >>logfile,' 17 17 VLA:E3 25.0 m -107.37.02.8 +33.54.00.5' 00270 print >>logfile,'' 00271 print >>logfile,'********** Regression ***********' 00272 print >>logfile,'* *' 00273 if (diff_mod<0.05): print >>logfile,'* Passed Model data test' 00274 print >>logfile,'* Model data mean'+str(thistest_mod)+','+str(model_amp) 00275 00276 if (diff_mod<0.05): 00277 regstate=True 00278 print >>logfile,'---' 00279 print >>logfile,'Passed Regression test for B0319' 00280 print >>logfile,'---' 00281 else: 00282 regstate=False 00283 print >>logfile,'----FAILED Regression test for B0319' 00284 print >>logfile,'*********************************' 00285 00286 print >>logfile,'' 00287 print >>logfile,'' 00288 print >>logfile,'********* Benchmarking *****************' 00289 print >>logfile,'* *' 00290 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00291 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00292 #print >>logfile,'Processing rate MB/s was: '+str(5.0/(endTime - startTime) 00293 print >>logfile,'* Breakdown: *' 00294 print >>logfile,'* import time was: '+str(importtime-startTime) 00295 print >>logfile,'* split time was: '+str(splittime-importtime) 00296 print >>logfile,'* plotant time was: '+str(plotanttime-splittime) 00297 print >>logfile,'* plotuvdist time was: '+str(plotuvdisttime-plotanttime) 00298 print >>logfile,'* plottime time was: '+str(plottimetime-plotuvdisttime) 00299 print >>logfile,'* plotchannel time was: '+str(plotchanneltime-plottimetime) 00300 print >>logfile,'* plotflag time was: '+str(plotflagtime-plotchanneltime) 00301 print >>logfile,'* plotaverage time was: '+str(plotavertime-plotflagtime) 00302 print >>logfile,'* calibrate time was: '+str(calibratetime-plotavertime) 00303 print >>logfile,'* plotcompare time was: '+str(plotcomparetime-calibratetime) 00304 print >>logfile,'* plotcal time was: '+str(plotcaltime-plotcomparetime) 00305 print >>logfile,'* plotcorr time was: '+str(plotcorrectedtime-plotcaltime) 00306 print >>logfile,'* uvmodelfit time was: '+str(uvmodelfittime-plotcorrectedtime) 00307 print >>logfile,'* ft time was: '+str(fttime-uvmodelfittime) 00308 print >>logfile,'* plotmodel time was: '+str(plotmodeltime-fttime) 00309 print >>logfile,'*****************************************' 00310 00311 logfile.close()
1.8.0
