|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ########################################################################## 00002 # # 00003 # Script for LeoRing EVLA Widar-0 Demo Science 20cm Data # 00004 # This is a 1 sub-band 256 channel OSRO mode dataset # 00005 # # 00006 # Created STM 2010-01-29 (Release R0) Regressionized Demo Script # 00007 # # 00008 ########################################################################## 00009 import time 00010 import os 00011 import pickle 00012 00013 #### Function needed to invert phase of visibilties 00014 # do scan-by-scan, spw-by-spw (ddids) visibility phase flip on requested column 00015 # assumes pylap is pl, and table component available 00016 # For column, specify 'DATA', 'CORRECTED_DATA' or 'MODEL_DATA' 00017 # Usage: signflip('data.ms','CORRECTED_DATA') 00018 def signflip (vis=None,column=None): 00019 00020 tb.open(vis,nomodify=False) 00021 00022 scancol=tb.getcol('SCAN_NUMBER') 00023 minscan=min(scancol) 00024 maxscan=max(scancol) 00025 00026 ddidcol=tb.getcol('DATA_DESC_ID') 00027 minddid=min(ddidcol) 00028 maxddid=max(ddidcol) 00029 00030 for scan in range(minscan,maxscan+1): 00031 for ddid in range(minddid,maxddid+1): 00032 st=tb.query('SCAN_NUMBER=='+str(scan)+' && DATA_DESC_ID=='+str(ddid)); 00033 if (st.nrows()>0): 00034 print 'scan=',scan, ' ddid=',ddid, ' rows=',st.nrows() 00035 d=tb.getcol(column) 00036 d=pl.conjugate(d) 00037 tb.putcol(column,d) 00038 st.done() 00039 tb.close() 00040 00041 ########################################################################## 00042 # # 00043 # Cleanup then run # 00044 # # 00045 ########################################################################## 00046 #### Clear previous runs 00047 os.system('rm -rf leo2pt_regression.*') 00048 00049 #------------------------------------------------------------------------- 00050 #Some variables 00051 #------------------------------------------------------------------------- 00052 regressname = 'LeoRing' 00053 prefix = 'leo2pt_regression' 00054 scriptprefix = 'run_leo2pt_regression' 00055 scriptvers = '20100129' 00056 00057 #------------------------------------------------------------------------- 00058 #Open an output logfile 00059 #------------------------------------------------------------------------- 00060 import datetime 00061 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00062 00063 outfile='out.'+prefix+'.'+datestring+'.log' 00064 logfile=open(outfile,'w') 00065 print 'Opened log file '+outfile 00066 00067 print >>logfile,'Processing SDM data for project '+regressname 00068 print >>logfile,'Script '+scriptvers 00069 00070 # Some date and version info 00071 myvers = casalog.version() 00072 myuser = os.getenv('USER') 00073 myhost = str( os.getenv('HOST') ) 00074 mycwd = os.getcwd() 00075 myos = os.uname() 00076 00077 # Print version info to outfile 00078 print >>logfile,'CASA '+myvers+' on host '+myhost+' running '+str(os.uname()[0]+os.uname()[2]) 00079 print >>logfile,'at '+datestring 00080 print >>logfile,'' 00081 print >>logfile,'' 00082 00083 stagename = [] 00084 stagetime = [] 00085 00086 startTime=time.time() 00087 startProc=time.clock() 00088 00089 prevTime = startTime 00090 00091 #------------------------------------------------------------------------- 00092 # Import, no autocorrs, then flag unaveraged data 00093 #------------------------------------------------------------------------- 00094 # Note: this raw (1s) SDM is available for download from 00095 # ftp://ftp.aoc.nrao.edu/pub/smyers/LeoRing/leo2pt.55183.452640752315.tgz 00096 # or 00097 # http://casa.nrao.edu/Data/EVLA/LeoRing/leo2pt.55183.452640752315.tgz 00098 print '---ImportASDM---' 00099 sdmfile = 'leo2pt.55183.452640752315' 00100 msfile = 'leo2pt_regression.ms' 00101 importasdm(asdm='leo2pt.55183.452640752315',vis='leo2pt_regression.ms',ocorr_mode='co') 00102 00103 # Some timing and file size stuff 00104 currTime=time.time() 00105 stagetime.append(currTime-prevTime) 00106 stagename.append('import') 00107 fillTime = currTime-prevTime 00108 prevTime = currTime 00109 # get dataset sizes in MB (1024x1024 bytes!) 00110 f = os.popen('du -ms '+sdmfile) 00111 fstr = f.readline() 00112 f.close() 00113 sdmsize = float( fstr.split("\t")[0] ) 00114 f = os.popen('du -ms '+msfile) 00115 fstr = f.readline() 00116 f.close() 00117 fillsize = float( fstr.split("\t")[0] ) 00118 00119 # Now flag on the unaveraged data 00120 # 00121 print '---Flagdata (Off-source timeranges)---' 00122 cmd = ["field='0' timerange='<10:52:02'", 00123 "field='4' timerange='13:44:38~13:47:57'"] 00124 tflagdata(vis='leo2pt_regression.ms',mode='list', inpfile=cmd) 00125 #flagdata(vis='leo2pt_regression.ms',mode='manualflag',field=['0','4'],timerange=['<10:52:02','13:44:38~13:47:57']) 00126 currTime=time.time() 00127 stagetime.append(currTime-prevTime) 00128 stagename.append('flagtime') 00129 prevTime = currTime 00130 00131 # Antenna 5, 15 supposedly bad 00132 print '---Flagdata (bad antennas)---' 00133 cmd = ["antenna='ea05'","antenna='ea15'"] 00134 tflagdata(vis='leo2pt_regression.ms',mode='list', inpfile=cmd) 00135 #flagdata(vis='leo2pt_regression.ms',mode='manualflag',antenna=['ea05','ea15']) 00136 currTime=time.time() 00137 stagetime.append(currTime-prevTime) 00138 stagename.append('flagants') 00139 prevTime = currTime 00140 00141 # 00142 print '---Flagdata (quack)---' 00143 tflagdata(vis='leo2pt_regression.ms',mode='quack',quackinterval=20.0) 00144 currTime=time.time() 00145 stagetime.append(currTime-prevTime) 00146 stagename.append('quack') 00147 prevTime = currTime 00148 00149 # 00150 print '---Flagdata (clip zeroes)---' 00151 tflagdata(vis='leo2pt_regression.ms',mode='clip', clipzeros=True, clipoutside=False, correlation='ABS_RR,LL') 00152 #flagdata(vis='leo2pt_regression.ms',mode='manualflag',clipminmax=[0.0,1.0E-8],clipoutside=F,clipexpr='ABS RR') 00153 #flagdata(vis='leo2pt_regression.ms',mode='manualflag',clipminmax=[0.0,1.0E-8],clipoutside=F,clipexpr='ABS LL') 00154 currTime=time.time() 00155 stagetime.append(currTime-prevTime) 00156 stagename.append('clipdata') 00157 prevTime = currTime 00158 00159 # Now split the data averaging to 10s 00160 # 00161 print '---Split (10s average)---' 00162 datafile = 'leo2pt_regression.split10sec.ms' 00163 split(vis='leo2pt_regression.ms',outputvis='leo2pt_regression.split10sec.ms',datacolumn='data',timebin='10s') 00164 currTime=time.time() 00165 stagetime.append(currTime-prevTime) 00166 stagename.append('split') 00167 splitTime = currTime-prevTime 00168 prevTime = currTime 00169 # get dataset sizes in MB (1024x1024 bytes!) 00170 f = os.popen('du -ms '+datafile) 00171 fstr = f.readline() 00172 f.close() 00173 splitsize = float( fstr.split("\t")[0] ) 00174 00175 # Change from EVLA to VLA 00176 print '---Vishead---' 00177 vishead(vis='leo2pt_regression.split10sec.ms',mode='put',hdkey='telescope',hdvalue='VLA') 00178 00179 # Create and initialize scratch columns 00180 print '---Clearcal---' 00181 clearcal(vis='leo2pt_regression.split10sec.ms') 00182 currTime=time.time() 00183 stagetime.append(currTime-prevTime) 00184 stagename.append('clearcal') 00185 clearTime = currTime-prevTime 00186 prevTime = currTime 00187 # get dataset sizes in MB (1024x1024 bytes!) 00188 f = os.popen('du -ms '+datafile) 00189 fstr = f.readline() 00190 f.close() 00191 clearsize = float( fstr.split("\t")[0] ) 00192 00193 print '---Listobs---' 00194 listobs('leo2pt_regression.split10sec.ms') 00195 currTime=time.time() 00196 stagetime.append(currTime-prevTime) 00197 stagename.append('listobs') 00198 prevTime = currTime 00199 00200 #### first gaincal on bandpass calibrator (temp) 00201 print '---Gaincal 1 (phase-only)---' 00202 gaincal(vis='leo2pt_regression.split10sec.ms',caltable='leo2pt_regression.gcal1', 00203 field='4',gaintype='G',calmode='p',refant='ea02', 00204 solint='int',spw='0:118~137') 00205 currTime=time.time() 00206 stagetime.append(currTime-prevTime) 00207 stagename.append('gaincal1') 00208 prevTime = currTime 00209 00210 plotcal(caltable='leo2pt_regression.gcal1', yaxis='phase', 00211 iteration='antenna',subplot=321,antenna='0~5', 00212 showgui=F,figfile='leo2pt_regression.gcal1.p1.png') 00213 00214 plotcal(caltable='leo2pt_regression.gcal1', yaxis='phase', 00215 iteration='antenna',subplot=321,antenna='6~11', 00216 showgui=F,figfile='leo2pt_regression.gcal1.p2.png') 00217 currTime=time.time() 00218 stagetime.append(currTime-prevTime) 00219 stagename.append('plot gcal1') 00220 prevTime = currTime 00221 00222 #### bandpass 00223 print '---Bandpass---' 00224 bandpass(vis='leo2pt_regression.split10sec.ms',caltable='leo2pt_regression.bcal', 00225 field='4',spw='',gaintable='leo2pt_regression.gcal1', 00226 interp=['nearest'], 00227 bandtype='B',solint='inf',combine='scan',refant='ea02',solnorm=True) 00228 currTime=time.time() 00229 stagetime.append(currTime-prevTime) 00230 stagename.append('bandpass1') 00231 prevTime = currTime 00232 00233 plotcal(caltable='leo2pt_regression.bcal',xaxis='freq',yaxis='amp',spw='', 00234 iteration='antenna',subplot=321,antenna='0~5', 00235 showgui=F,figfile='leo2pt_regression.bcal.p1.png') 00236 00237 plotcal(caltable='leo2pt_regression.bcal',xaxis='freq',yaxis='amp',spw='', 00238 iteration='antenna',subplot=321,antenna='6~11', 00239 showgui=F,figfile='leo2pt_regression.bcal.p2.png') 00240 currTime=time.time() 00241 stagetime.append(currTime-prevTime) 00242 stagename.append('plot bcal1') 00243 prevTime = currTime 00244 00245 #### final phasecal 00246 print '---Gaincal 2 (phase-only)---' 00247 gaincal(vis='leo2pt_regression.split10sec.ms',caltable='leo2pt_regression.gcal2', 00248 gaintable='leo2pt_regression.bcal', interp='nearest', 00249 field='1,4',gaintype='G',calmode='p',refant='ea02', 00250 solint='int',spw='0:30~224') 00251 currTime=time.time() 00252 stagetime.append(currTime-prevTime) 00253 stagename.append('gaincal2') 00254 prevTime = currTime 00255 00256 plotcal(caltable='leo2pt_regression.gcal2', yaxis='phase', 00257 iteration='antenna',subplot=321,antenna='0~5', 00258 showgui=F,figfile='leo2pt_regression.gcal2.p1.png') 00259 00260 plotcal(caltable='leo2pt_regression.gcal2', yaxis='phase', 00261 iteration='antenna',subplot=321,antenna='6~11', 00262 showgui=F,figfile='leo2pt_regression.gcal2.p2.png') 00263 currTime=time.time() 00264 stagetime.append(currTime-prevTime) 00265 stagename.append('plot gcal2') 00266 prevTime = currTime 00267 00268 #### final ampcal 00269 00270 print '---Table (change field name)---' 00271 # For this to work had to change the Field name for field 4 00272 # from J1331+3030 to 1331+305 00273 tb.open('leo2pt_regression.split10sec.ms/FIELD',nomodify=False) 00274 st=tb.selectrows(4) 00275 st.putcol('NAME','1331+305') 00276 st.done() 00277 tb.close() 00278 00279 # Calculates 14.7461Jy 00280 print '---Setjy---' 00281 setjy(vis='leo2pt_regression.split10sec.ms',field='4',modimage='3C286_L.im',scalebychan=False,standard='Perley-Taylor 99') 00282 currTime=time.time() 00283 stagetime.append(currTime-prevTime) 00284 stagename.append('setjy') 00285 prevTime = currTime 00286 00287 print '---Gaincal 3 (amp-only)---' 00288 gaincal(vis='leo2pt_regression.split10sec.ms',caltable='leo2pt_regression.gcal3', 00289 gaintable=['leo2pt_regression.bcal','leo2pt_regression.gcal2'], 00290 interp=['nearest','nearest'], 00291 field='1,4',gaintype='G',calmode='a',refant='ea02', 00292 solint='inf',combine='',spw='0:30~224') 00293 currTime=time.time() 00294 stagetime.append(currTime-prevTime) 00295 stagename.append('gaincal3') 00296 prevTime = currTime 00297 00298 plotcal(caltable='leo2pt_regression.gcal3', yaxis='amp', 00299 iteration='antenna',subplot=321,antenna='0~5', 00300 showgui=F,figfile='leo2pt_regression.gcal3.p1.png') 00301 00302 plotcal(caltable='leo2pt_regression.gcal3', yaxis='amp', 00303 iteration='antenna',subplot=321,antenna='6~11', 00304 showgui=F,figfile='leo2pt_regression.gcal3.p2.png') 00305 currTime=time.time() 00306 stagetime.append(currTime-prevTime) 00307 stagename.append('plot gcal3') 00308 prevTime = currTime 00309 00310 print '---Fluxscale---' 00311 fluxscale(vis='leo2pt_regression.split10sec.ms',caltable='leo2pt_regression.gcal3', 00312 fluxtable='leo2pt_regression.fscale',reference='4') 00313 currTime=time.time() 00314 stagetime.append(currTime-prevTime) 00315 stagename.append('fluxscale') 00316 prevTime = currTime 00317 00318 plotcal(caltable='leo2pt_regression.fscale', yaxis='amp', 00319 iteration='antenna',subplot=321,antenna='0~5', 00320 showgui=F,figfile='leo2pt_regression.fscale.p1.png') 00321 00322 plotcal(caltable='leo2pt_regression.fscale', yaxis='amp', 00323 iteration='antenna',subplot=321,antenna='6~11', 00324 showgui=F,figfile='leo2pt_regression.fscale.p2.png') 00325 currTime=time.time() 00326 stagetime.append(currTime-prevTime) 00327 stagename.append('plot fluxscale') 00328 prevTime = currTime 00329 00330 #### final applycal 00331 print '---Applycal (field 1)---' 00332 applycal(vis='leo2pt_regression.split10sec.ms',spw='', field='1', 00333 gaintable=['leo2pt_regression.bcal','leo2pt_regression.gcal2','leo2pt_regression.fscale'], 00334 interp=['nearest','nearest','nearest'],gainfield=['4','1','1']) 00335 currTime=time.time() 00336 stagetime.append(currTime-prevTime) 00337 stagename.append('applycal 1') 00338 prevTime = currTime 00339 00340 print '---Applycal (field 2,3)---' 00341 # transfer amp+phase from field 1 to field 2,3 00342 applycal(vis='leo2pt_regression.split10sec.ms',spw='', field='2,3', 00343 gaintable=['leo2pt_regression.bcal','leo2pt_regression.gcal2','leo2pt_regression.fscale'], 00344 interp=['nearest','linear','linear'],gainfield=['4','1','1']) 00345 currTime=time.time() 00346 stagetime.append(currTime-prevTime) 00347 stagename.append('applycal 2,3') 00348 prevTime = currTime 00349 00350 print '---Applycal (field 4)---' 00351 applycal(vis='leo2pt_regression.split10sec.ms',spw='', field='4', 00352 gaintable=['leo2pt_regression.bcal','leo2pt_regression.gcal2','leo2pt_regression.fscale'], 00353 interp=['nearest','nearest','nearest'],gainfield=['4','4','4']) 00354 currTime=time.time() 00355 stagetime.append(currTime-prevTime) 00356 stagename.append('applycal 4') 00357 prevTime = currTime 00358 00359 #### split calibrated data 00360 print '---Split (field 1)---' 00361 split(vis='leo2pt_regression.split10sec.ms',outputvis='leo2pt_regression.field1.ms', 00362 datacolumn='corrected',field='1') 00363 currTime=time.time() 00364 stagetime.append(currTime-prevTime) 00365 stagename.append('split field 1') 00366 prevTime = currTime 00367 00368 print '---Split (field 2)---' 00369 split(vis='leo2pt_regression.split10sec.ms',outputvis='leo2pt_regression.field2.ms', 00370 datacolumn='corrected',field='2') 00371 currTime=time.time() 00372 stagetime.append(currTime-prevTime) 00373 stagename.append('split field 2') 00374 prevTime = currTime 00375 00376 print '---Split (field 3)---' 00377 split(vis='leo2pt_regression.split10sec.ms',outputvis='leo2pt_regression.field3.ms', 00378 datacolumn='corrected',field='3') 00379 currTime=time.time() 00380 stagetime.append(currTime-prevTime) 00381 stagename.append('split field 3') 00382 prevTime = currTime 00383 00384 print '---Split (field 4)---' 00385 split(vis='leo2pt_regression.split10sec.ms',outputvis='leo2pt_regression.field4.ms', 00386 datacolumn='corrected',field='4') 00387 currTime=time.time() 00388 stagetime.append(currTime-prevTime) 00389 stagename.append('split field 4') 00390 prevTime = currTime 00391 00392 print '---Split (all)---' 00393 split(vis='leo2pt_regression.split10sec.ms',outputvis='leo2pt_regression.all.ms', 00394 datacolumn='corrected',field='1~4') 00395 currTime=time.time() 00396 stagetime.append(currTime-prevTime) 00397 stagename.append('split all') 00398 prevTime = currTime 00399 00400 #### flip signs 00401 print '---Signflip (field 1)---' 00402 signflip('leo2pt_regression.field1.ms','DATA') 00403 currTime=time.time() 00404 stagetime.append(currTime-prevTime) 00405 stagename.append('signflip 1') 00406 prevTime = currTime 00407 00408 print '---Signflip (field 2)---' 00409 signflip('leo2pt_regression.field2.ms','DATA') 00410 currTime=time.time() 00411 stagetime.append(currTime-prevTime) 00412 stagename.append('signflip 2') 00413 prevTime = currTime 00414 00415 print '---Signflip (field 3)---' 00416 signflip('leo2pt_regression.field3.ms','DATA') 00417 currTime=time.time() 00418 stagetime.append(currTime-prevTime) 00419 stagename.append('signflip 3') 00420 prevTime = currTime 00421 00422 print '---Signflip (field 4)---' 00423 signflip('leo2pt_regression.field4.ms','DATA') 00424 currTime=time.time() 00425 stagetime.append(currTime-prevTime) 00426 stagename.append('signflip 4') 00427 prevTime = currTime 00428 00429 print '---Signflip (all)---' 00430 signflip('leo2pt_regression.all.ms','DATA') 00431 currTime=time.time() 00432 stagetime.append(currTime-prevTime) 00433 stagename.append('signflip all') 00434 prevTime = currTime 00435 00436 print '---ExportUVFITS---' 00437 fitsfile = 'leo2pt_regression.all.uvfits' 00438 exportuvfits(vis='leo2pt_regression.all.ms',fitsfile='leo2pt_regression.all.uvfits', 00439 datacolumn='data',multisource=True,async=False) 00440 currTime=time.time() 00441 stagetime.append(currTime-prevTime) 00442 stagename.append('exportuvfits') 00443 exportTime = currTime-prevTime 00444 prevTime = currTime 00445 # get dataset sizes in MB (1024x1024 bytes!) 00446 f = os.popen('du -ms leo2pt_regression.all.ms') 00447 fstr = f.readline() 00448 f.close() 00449 msallsize = float( fstr.split("\t")[0] ) 00450 # get dataset sizes in MB (1024x1024 bytes!) 00451 f = os.popen('du -ms '+fitsfile) 00452 fstr = f.readline() 00453 f.close() 00454 fitssize = float( fstr.split("\t")[0] ) 00455 # 00456 ########################################################################## 00457 # 00458 # Get MS stats 00459 # 00460 print '---Visstat (J1042+1203)---' 00461 visstat_cal=visstat('leo2pt_regression.field1.ms') 00462 vismean_cal=visstat_cal['DATA']['mean'] 00463 print "Found vis mean on J1042+1203 = "+str(vismean_cal) 00464 currTime=time.time() 00465 stagetime.append(currTime-prevTime) 00466 stagename.append('visstat fld 1') 00467 prevTime = currTime 00468 00469 print '---Visstat (1331+305)---' 00470 visstat_3c286=visstat('leo2pt_regression.field4.ms') 00471 vismean_3c286=visstat_3c286['DATA']['mean'] 00472 print "Found vis mean on 1331+305 = "+str(vismean_3c286) 00473 currTime=time.time() 00474 stagetime.append(currTime-prevTime) 00475 stagename.append('visstat fld 4') 00476 prevTime = currTime 00477 # 00478 ########################################################################## 00479 #First cvel to bary 00480 00481 print '---Cvel (field 2)---' 00482 cvel(vis='leo2pt_regression.field2.ms', 00483 outputvis ='leo2pt_regression.field2_chanbary.ms', 00484 mode='channel',nchan=-1,start=0,width=1, 00485 interpolation='nearest', 00486 phasecenter='', 00487 spw='', 00488 restfreq='1420405751.786Hz', 00489 outframe='BARY') 00490 currTime=time.time() 00491 stagetime.append(currTime-prevTime) 00492 stagename.append('cvel fld 2') 00493 prevTime = currTime 00494 00495 print '---Cvel (field 3)---' 00496 cvel(vis='leo2pt_regression.field3.ms', 00497 outputvis ='leo2pt_regression.field3_chanbary.ms', 00498 mode='channel',nchan=-1,start=0,width=1, 00499 interpolation='nearest', 00500 phasecenter='', 00501 spw='', 00502 restfreq='1420405751.786Hz', 00503 outframe='BARY') 00504 currTime=time.time() 00505 stagetime.append(currTime-prevTime) 00506 stagename.append('cvel fld 3') 00507 prevTime = currTime 00508 00509 #### Combine Leo-1 and Leo-2 fields 00510 # 2 NONE Leo-1 10:47:22.0000 +12.16.38.0000 J2000 2 301290 00511 # 3 NONE Leo-2 10:46:45.0000 +11.50.38.0000 J2000 3 300960 00512 print '---Concat field 2 and 3 (Leo-1 and Leo-2)---' 00513 concat(vis = ['leo2pt_regression.field2_chanbary.ms','leo2pt_regression.field3_chanbary.ms'], 00514 freqtol='150Hz', 00515 concatvis = 'leo2pt_regression.concat_chanbary.ms') 00516 currTime=time.time() 00517 stagetime.append(currTime-prevTime) 00518 stagename.append('concat') 00519 prevTime = currTime 00520 00521 #listobs('leo2pt_regression.concat_chanbary.ms') 00522 00523 # combine with center at 10:47:03.50 +12.03.38.0000 00524 print '---Clean (dirty mosaic)---' 00525 clean(vis='leo2pt_regression.concat_chanbary.ms', 00526 imagename='leo2pt_regression.concat_chanbary.dirtymos', 00527 imsize=[500,500],cell=[12.0,12.0],niter=0, 00528 mode='channel',nchan=-1,start=0,width=1, 00529 interpolation='nearest', 00530 psfmode='hogbom',imagermode='mosaic',ftmachine='ft', 00531 minpb=0.1,pbcor=False, 00532 phasecenter='J2000 10h47m03.50s 12d03m38.0s', 00533 restfreq='1420405751.786Hz',weighting='natural') 00534 currTime=time.time() 00535 stagetime.append(currTime-prevTime) 00536 stagename.append('dirty') 00537 prevTime = currTime 00538 00539 #viewer('leo2pt_regression.concat_chanbary.dirtymos.image') 00540 # 00541 # Statistics on dirty image cube 00542 # 00543 print '--ImStat (Dirty cube)--' 00544 dirtystat = imstat('leo2pt_regression.concat_chanbary.dirtymos.image') 00545 print "Found dirty image max = "+str(dirtystat['max'][0]) 00546 print "Found dirty image rms = "+str(dirtystat['sigma'][0]) 00547 00548 dirtyoff = imstat('leo2pt_regression.concat_chanbary.dirtymos.image', 00549 chans='',box='175,190,210,230') 00550 print "Found off-source dirty image rms = "+str(dirtyoff['sigma'][0]) 00551 00552 currTime=time.time() 00553 stagetime.append(currTime-prevTime) 00554 stagename.append('imstat dirty') 00555 prevTime = currTime 00556 00557 #------------------------------------------- 00558 # Clean all at once 00559 print '---Clean (ft mosaic)---' 00560 clean(vis='leo2pt_regression.concat_chanbary.ms', 00561 imagename='leo2pt_regression.concat_chanbary.cleanmosft', 00562 imsize=[500,500],cell=[12.0,12.0], 00563 niter=40000,threshold='20mJy', 00564 mode='channel',nchan=256,start=0,width=1, 00565 interpolation='nearest', 00566 psfmode='hogbom',imagermode='mosaic',ftmachine='ft', 00567 minpb=0.1,pbcor=False, 00568 phasecenter='J2000 10h47m03.50s 12d03m38.0s', 00569 restfreq='1420405751.786Hz',weighting='natural') 00570 currTime=time.time() 00571 stagetime.append(currTime-prevTime) 00572 stagename.append('clean') 00573 prevTime = currTime 00574 00575 #Fitted beam used in restoration: 86.4579 by 48.7073 (arcsec) at pa 178.107 (deg) 00576 00577 #viewer('leo2pt_regression.concat_chanbary.cleanmosft.image') 00578 # 00579 # Statistics on clean image cube 00580 # 00581 print '--ImStat (Clean cube)--' 00582 cleanstat = imstat('leo2pt_regression.concat_chanbary.cleanmosft.image') 00583 print "Found clean image max = "+str(cleanstat['max'][0]) 00584 print "Found clean image rms = "+str(cleanstat['sigma'][0]) 00585 00586 offstat = imstat('leo2pt_regression.concat_chanbary.cleanmosft.image', 00587 chans='',box='175,190,210,230') 00588 print "Found off-source clean image rms = "+str(offstat['sigma'][0]) 00589 00590 # Statistics on clean model 00591 print '--ImStat (Clean model)--' 00592 modstat = imstat('leo2pt_regression.concat_chanbary.cleanmosft.model') 00593 print "Found total model flux = "+str(modstat['sum'][0]) 00594 00595 currTime=time.time() 00596 stagetime.append(currTime-prevTime) 00597 stagename.append('imstat clean') 00598 prevTime = currTime 00599 00600 #------------------------------------------- 00601 # Clean all at once with uvtaper 00602 #3klambda: 00603 #Fitted beam used in restoration: 137.297 by 58.0881 (arcsec) at pa 158.811 (deg) 00604 00605 #2.2klambda 00606 #Fitted beam used in restoration: 153.274 by 66.2216 (arcsec) at pa 160.137 (deg) 00607 00608 #1.6klambda 00609 #Fitted beam used in restoration: 180.444 by 79.4114 (arcsec) at pa 156.049 (deg) 00610 00611 print '---Clean (ft mosaic tapered)---' 00612 clean(vis='leo2pt_regression.concat_chanbary.ms', 00613 imagename='leo2pt_regression.concat_chanbary.taper1600.cleanmosft', 00614 imsize=[400,400],cell=[20.0,20.0], 00615 niter=40000,threshold='25mJy', 00616 mode='channel',nchan=256,start=0,width=1, 00617 interpolation='nearest', 00618 psfmode='hogbom',imagermode='mosaic',ftmachine='ft', 00619 minpb=0.2,pbcor=False,cyclefactor=1.0, 00620 phasecenter='J2000 10h47m03.50s 12d03m38.0s', 00621 uvtaper=T,outertaper=['1600lambda'], 00622 restfreq='1420405751.786Hz',weighting='natural') 00623 currTime=time.time() 00624 stagetime.append(currTime-prevTime) 00625 stagename.append('clean tapered') 00626 prevTime = currTime 00627 00628 #viewer('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image') 00629 # 00630 # Statistics on clean image cube 00631 # 00632 print '--ImStat (Clean tapered image)--' 00633 taperstat = imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image') 00634 print "Found tapered image max = "+str(taperstat['max'][0]) 00635 print "Found tapered image rms = "+str(taperstat['sigma'][0]) 00636 00637 taperoff = imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image', 00638 chans='',box='175,190,210,230') 00639 print "Found off-source tapered image rms = "+str(taperoff['sigma'][0]) 00640 00641 # Statistics on tapered clean model 00642 print '--ImStat (Clean tapered model)--' 00643 tapermod = imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.model') 00644 print "Found total tapered model flux = "+str(tapermod['sum'][0]) 00645 00646 currTime=time.time() 00647 stagetime.append(currTime-prevTime) 00648 stagename.append('imstat tapered') 00649 prevTime = currTime 00650 00651 ########################################################################## 00652 #Moment 0 00653 print '---Immoments (0)---' 00654 immoments(imagename='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image', 00655 moments=[0],axis='spectral', 00656 excludepix=[-1000.,0.015], 00657 outfile='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.moment0') 00658 currTime=time.time() 00659 stagetime.append(currTime-prevTime) 00660 stagename.append('moment0') 00661 prevTime = currTime 00662 00663 # Remove the pixel mask (for better display) 00664 ia.open('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.moment0') 00665 ia.maskhandler('delete','mask0') 00666 ia.close() 00667 00668 #Moment 1 00669 print '---Immoments (1)---' 00670 immoments(imagename='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image', 00671 moments=[1],axis='spectral', 00672 excludepix=[-1000.,0.025], 00673 outfile='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.moment1') 00674 currTime=time.time() 00675 stagetime.append(currTime-prevTime) 00676 stagename.append('moment1') 00677 prevTime = currTime 00678 00679 #viewer('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.moment0') 00680 00681 # 00682 # Statistics on moment images 00683 # 00684 print '--ImStat (Moment images)--' 00685 momzerostat=imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.moment0') 00686 try: 00687 print "Found moment 0 max = "+str(momzerostat['max'][0]) 00688 print "Found moment 0 rms = "+str(momzerostat['rms'][0]) 00689 except: 00690 pass 00691 00692 momonestat=imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.moment1') 00693 try: 00694 print "Found moment 1 median = "+str(momonestat['median'][0]) 00695 except: 00696 pass 00697 00698 # Get velocities of first and last planes of cube 00699 print '---Immoments (chans 0 and 255)---' 00700 # Do a moment one on channel 0 to check that the indexing is right 00701 try: 00702 immoments(imagename='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image', 00703 moments=[1],includepix=[], 00704 chans='0', 00705 outfile='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.plane0.mom1') 00706 except: 00707 pass 00708 00709 # Do a moment one on channel 255 to check that the indexing is right 00710 try: 00711 immoments(imagename='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image', 00712 moments=[1],includepix=[], 00713 chans='255', 00714 outfile='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.plane255.mom1') 00715 except: 00716 pass 00717 00718 ia.open('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image') 00719 csys=ia.coordsys() 00720 vel0=0.0 00721 vel255=0.0 00722 00723 try: 00724 momoneplane0=imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.plane0.mom1') 00725 print "Found plane 0 moment 1 value = "+str(momoneplane0['median'][0]) 00726 except: 00727 pass 00728 00729 00730 try: 00731 momoneplane255=imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.plane255.mom1') 00732 print "Found plane 255 moment 1 value = "+str(momoneplane255['median'][0]) 00733 except: 00734 pass 00735 00736 if(type(momoneplane0)==bool): 00737 vel0=csys.frequencytovelocity(ia.toworld([0,0,0,0])['numeric'][3]) 00738 if(type(momoneplane255)==bool): 00739 vel255=csys.frequencytovelocity(ia.toworld([0,0,0,255])['numeric'][3]) 00740 00741 currTime=time.time() 00742 stagetime.append(currTime-prevTime) 00743 stagename.append('imstat moments') 00744 prevTime = currTime 00745 # 00746 ########################################################################## 00747 # 00748 # Manually correct for mosaic response pattern using .image/.flux images 00749 print '--ImMath (PBcor)--' 00750 immath(outfile='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.pbcor', 00751 imagename=['leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image', 00752 'leo2pt_regression.concat_chanbary.taper1600.cleanmosft.flux'], 00753 mode='evalexpr',expr='IM0/IM1') 00754 00755 # now pbcor the model, be careful to mask zeros 00756 immath(outfile='leo2pt_regression.concat_chanbary.taper1600.cleanmosft.pbcormod', 00757 imagename=['leo2pt_regression.concat_chanbary.taper1600.cleanmosft.model', 00758 'leo2pt_regression.concat_chanbary.taper1600.cleanmosft.flux'], 00759 mode='evalexpr',expr='IM0/IM1[IM1!=0.0]') 00760 currTime=time.time() 00761 stagetime.append(currTime-prevTime) 00762 stagename.append('immath pbcor') 00763 prevTime = currTime 00764 00765 # 00766 # Statistics on PBcor image cube 00767 # 00768 print '--ImStat (PBcor cube)--' 00769 pbcorstat = imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.pbcor') 00770 print "Found pbcor image max = "+str(pbcorstat['max'][0]) 00771 00772 pbcoroffstat = imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.pbcor', 00773 chans='',box='175,190,210,230') 00774 print "Found pbcor off-source rms = "+str(pbcoroffstat['sigma'][0]) 00775 00776 # 00777 # Statistics on PBcor model cube 00778 # 00779 print '--ImStat (PBcor model)--' 00780 pbcormodstat = imstat('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.pbcormod') 00781 print "Found total model flux = "+str(pbcormodstat['sum'][0]) 00782 00783 currTime=time.time() 00784 stagetime.append(currTime-prevTime) 00785 stagename.append('imstat pbcor') 00786 prevTime = currTime 00787 00788 # 00789 ########################################################################## 00790 # 00791 #DS9 shape files for viewer 00792 #From ALFALFA 00793 #j2000; text 10:46:45.7 11:49:12 # text={M96} color=blue 00794 #j2000; text 10:47:20.1 12:23:15 # text={205505} color=green 00795 #j2000; text 10:46:41.3 12:19:37 # text={202027} color=cyan 00796 #From EVLA 00797 #j2000; text 10:46:05.7 12:08:03 # text={A} color=magenta 00798 #j2000; text 10:47:56.6 12:22:18 # text={B} color=yellow 00799 #j2000; text 10:48:17.2 12:21:48 # text={C} color=red 00800 #j2000; text 10:46:51.4 11:54:53 # text={D} color=blue 00801 00802 #pixels: 00803 #M96 210,156 00804 #205505 00805 #202027 209,247 00806 #A 00807 #B 00808 #C 146,255 00809 00810 #---------------------------- 00811 #Make spectrum at M96 00812 print '---Imval (M96)---' 00813 myval = imval('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image',box='210,156,210,156') 00814 M96_sparr = myval['data'] 00815 00816 ia.open('leo2pt_regression.concat_chanbary.taper1600.cleanmosft.image') 00817 chlist = range(256) 00818 fqlist = [] 00819 #find ra,dec 00820 mys = ia.toworld([210,156,0,0],'s') 00821 M96_ra = mys['string'][0] 00822 M96_dec = mys['string'][1] 00823 #find spectral coordinates 00824 for i in chlist: 00825 myw = ia.toworld([210,156,0,i],'n') 00826 fqlist.append(myw['numeric'][3]) 00827 00828 ia.close() 00829 M96_fqlist = fqlist 00830 00831 fqarr = pl.array(fqlist) 00832 #convert to GHz from Hz 00833 M96_fqarr = fqarr*1.0E-9 00834 00835 print 'LeoRing spectrum at M96 : ',M96_ra,' ',M96_dec 00836 currTime=time.time() 00837 stagetime.append(currTime-prevTime) 00838 stagename.append('imval m96') 00839 prevTime = currTime 00840 00841 # Now use matplotlib to plot this spectrum 00842 print '---Matplotlib M96---' 00843 pl.plot(M96_fqarr, M96_sparr, linewidth=1.0) 00844 pl.ylim(-0.04,1.4) 00845 pl.xlim(23.6,26.7) 00846 pl.xlabel('Barycentric Frequency (GHz)') 00847 pl.ylabel('Flux Density (Jy/beam)') 00848 pl.title('EVLA Leo Ring at M96 = 10:46:45.7 11d49m12s') 00849 # Save plot 00850 pl.savefig('leo2pt_regression.M96_spectrum.png',format='png') 00851 pl.savefig('leo2pt_regression.M96_spectrum.eps',dpi=600,format='eps') 00852 pl.savefig('leo2pt_regression.M96_spectrum.pdf',dpi=600,format='pdf') 00853 pl.close() 00854 00855 #Write out spectrum to ascii file: 00856 outfile = 'leo2pt_regression.M96_spectrum.txt' 00857 specfile=open(outfile,'w') 00858 00859 print >>specfile,'#EVLA Leo Ring at M96 = 10:46:45.7 11d49m12s' 00860 print >>specfile,'#256 points: BARYvelocity(km/s) FluxDensity(Jy/beam)' 00861 for i in chlist: 00862 print >>specfile,'%10.7f %12.6f' % (M96_fqarr[i], M96_sparr[i]) 00863 00864 specfile.close() 00865 currTime=time.time() 00866 stagetime.append(currTime-prevTime) 00867 stagename.append('matplotlib m96') 00868 prevTime = currTime 00869 00870 #================================================================================ 00871 # Done 00872 endProc=time.clock() 00873 endTime=time.time() 00874 00875 #================================================================================ 00876 # Set up the regression values 00877 # Defaults 00878 dirty_image_max = 1.0 00879 dirty_offsrc_rms = 1.0 00880 clean_image_max = 1.0 00881 clean_offsrc_rms = 1.0 00882 taper_image_max = 1.0 00883 taper_offsrc_rms = 1.0 00884 clean_momentzero_max = 1.0 00885 clean_momentzero_rms = 1.0 00886 clean_momentone_median = 0.0 00887 clean_momentone_planezero = 0.0 00888 clean_momentone_planelast = 0.0 00889 vis_mean_cal = 1.0 00890 vis_mean_3c286 = 1.0 00891 model_sum = 1.0 00892 model_clean_sum = 1.0 00893 model_taper_sum = 1.0 00894 model_pbcor_sum = 1.0 00895 00896 # STM 2010-01-29 Version 3.0.1 build 10067 00897 testdate = '2010-01-29 (STM)' 00898 testvers = 'CASA Version 3.0.1 Rev 10067' 00899 dirty_image_max = 0.14097676 00900 dirty_offsrc_rms = 0.0129723958478 00901 clean_image_max = 0.120128452778 00902 clean_offsrc_rms = 0.00723671344488 00903 taper_image_max = 0.177342906594 00904 taper_offsrc_rms = 0.00835683070787 00905 clean_momentzero_max = 28.124853 00906 clean_momentzero_rms = 0.650879 00907 clean_momentone_median = 852.958618 00908 clean_momentone_planezero = 1089.798584 00909 clean_momentone_planelast = 669.364807 00910 vis_mean_cal = 3.565428 00911 vis_mean_3c286 = 14.837324 00912 model_clean_sum = 51.4040303375 00913 model_taper_sum = 51.293914702 00914 model_pbcor_sum = 73.011602 00915 00916 canonical = {} 00917 canonical['exist'] = True 00918 00919 canonical['date'] = testdate 00920 canonical['version'] = testvers 00921 canonical['user'] = 'smyers' 00922 canonical['host'] = 'rishi' 00923 canonical['cwd'] = '/home/rishi2/smyers/3.0.1/LeoRing/' 00924 print "Using internal regression from "+canonical['version']+" on "+canonical['date'] 00925 00926 canonical_results = {} 00927 canonical_results['dirty_image_max'] = {} 00928 canonical_results['dirty_image_max']['value'] = dirty_image_max 00929 canonical_results['dirty_image_offsrc_max'] = {} 00930 canonical_results['dirty_image_offsrc_max']['value'] = dirty_offsrc_rms 00931 canonical_results['clean_image_max'] = {} 00932 canonical_results['clean_image_max']['value'] = clean_image_max 00933 canonical_results['clean_image_offsrc_max'] = {} 00934 canonical_results['clean_image_offsrc_max']['value'] = clean_offsrc_rms 00935 canonical_results['taper_image_max'] = {} 00936 canonical_results['taper_image_max']['value'] = taper_image_max 00937 canonical_results['taper_image_offsrc_max'] = {} 00938 canonical_results['taper_image_offsrc_max']['value'] = taper_offsrc_rms 00939 canonical_results['clean_momentzero_max'] = {} 00940 canonical_results['clean_momentzero_max']['value'] = clean_momentzero_max 00941 canonical_results['clean_momentzero_rms'] = {} 00942 canonical_results['clean_momentzero_rms']['value'] = clean_momentzero_rms 00943 canonical_results['clean_momentone_median'] = {} 00944 canonical_results['clean_momentone_median']['value'] = clean_momentone_median 00945 canonical_results['clean_momentone_planezero'] = {} 00946 canonical_results['clean_momentone_planezero']['value'] = clean_momentone_planezero 00947 canonical_results['clean_momentone_planelast'] = {} 00948 canonical_results['clean_momentone_planelast']['value'] = clean_momentone_planelast 00949 00950 canonical_results['vis_mean_cal'] = {} 00951 canonical_results['vis_mean_cal']['value'] = vis_mean_cal 00952 canonical_results['vis_mean_3c286'] = {} 00953 canonical_results['vis_mean_3c286']['value'] = vis_mean_3c286 00954 00955 canonical_results['model_clean_sum'] = {} 00956 canonical_results['model_clean_sum']['value'] = model_clean_sum 00957 canonical_results['model_taper_sum'] = {} 00958 canonical_results['model_taper_sum']['value'] = model_taper_sum 00959 canonical_results['model_pbcor_sum'] = {} 00960 canonical_results['model_pbcor_sum']['value'] = model_pbcor_sum 00961 00962 canonical['results'] = canonical_results 00963 00964 print "Canonical Regression (default) from "+canonical['date'] 00965 00966 # 00967 # Try and load previous results from regression file 00968 # 00969 regression = {} 00970 regressfile = scriptprefix + '.pickle' 00971 prev_results = {} 00972 00973 try: 00974 fr = open(regressfile,'r') 00975 except: 00976 print "No previous regression results file "+regressfile 00977 regression['exist'] = False 00978 else: 00979 u = pickle.Unpickler(fr) 00980 regression = u.load() 00981 fr.close() 00982 print "Regression results filled from "+regressfile 00983 print "Regression from version "+regression['version']+" on "+regression['date'] 00984 regression['exist'] = True 00985 00986 prev_results = regression['results'] 00987 00988 # 00989 ########################################################################## 00990 # Calculate test values 00991 ########################################################################## 00992 # 00993 print '--Calculate Results--' 00994 print '' 00995 00996 try: 00997 dirtymax = dirtystat['max'][0] 00998 except: 00999 dirtymax = 0.0 01000 01001 try: 01002 dirtyoffrms = dirtyoff['sigma'][0] 01003 except: 01004 dirtyoffrms = 0.0 01005 01006 try: 01007 cleanmax = cleanstat['max'][0] 01008 except: 01009 cleanmax = 0.0 01010 01011 try: 01012 cleanoffrms = offstat['sigma'][0] 01013 except: 01014 cleanoffrms = 0.0 01015 01016 try: 01017 tapermax = taperstat['max'][0] 01018 except: 01019 tapermax = 0.0 01020 01021 try: 01022 taperoffrms = taperoff['sigma'][0] 01023 except: 01024 taperoffrms = 0.0 01025 01026 try: 01027 momzero_max = momzerostat['max'][0] 01028 except: 01029 momzero_max = 0.0 01030 01031 try: 01032 momzero_rms = momzerostat['rms'][0] 01033 except: 01034 momzero_rms = 0.0 01035 01036 try: 01037 momone_median = momonestat['median'][0] 01038 except: 01039 momone_median = 0.0 01040 01041 try: 01042 momone_plane0 = momoneplane0['median'][0] 01043 except: 01044 momone_plane0 = vel0 01045 01046 try: 01047 momone_plane255 = momoneplane255['median'][0] 01048 except: 01049 momone_plane255 = vel255 01050 01051 try: 01052 cleanmodflux = modstat['sum'][0] 01053 except: 01054 cleanmodflux = 0.0 01055 01056 try: 01057 tapermodflux = tapermod['sum'][0] 01058 except: 01059 tapermodflux = 0.0 01060 01061 try: 01062 pbcormodflux = pbcormodstat['sum'][0] 01063 except: 01064 pbcormodflux = 0.0 01065 01066 # 01067 # Store results in dictionary 01068 # 01069 new_regression = {} 01070 01071 # Some date and version info 01072 import datetime 01073 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 01074 01075 # System info 01076 myvers = casalog.version() 01077 try: 01078 myuser = os.getlogin() 01079 except: 01080 myuser = os.getenv('USER') 01081 #myhost = str( os.getenv('HOST') ) 01082 myuname = os.uname() 01083 myhost = myuname[1] 01084 myos = myuname[0]+' '+myuname[2]+' '+myuname[4] 01085 mycwd = os.getcwd() 01086 mypath = os.environ.get('CASAPATH') 01087 01088 mydataname = 'LeoRing 18-Dec-2009 EVLA' 01089 mydataset = sdmfile 01090 01091 # Save info in regression dictionary 01092 new_regression['date'] = datestring 01093 new_regression['version'] = myvers 01094 new_regression['user'] = myuser 01095 new_regression['host'] = myhost 01096 new_regression['cwd'] = mycwd 01097 new_regression['os'] = myos 01098 new_regression['uname'] = myuname 01099 new_regression['aipspath'] = mypath 01100 01101 new_regression['dataname'] = mydataname 01102 new_regression['dataset'] = mydataset 01103 01104 # Fill results 01105 # Note that 'op' tells what to do for the diff : 01106 # 'divf' = abs( new - prev )/prev 01107 # 'diff' = new - prev 01108 01109 results = {} 01110 01111 op = 'divf' 01112 tol = 0.08 01113 results['dirty_image_max'] = {} 01114 results['dirty_image_max']['name'] = 'Dirty image max' 01115 results['dirty_image_max']['value'] = dirtymax 01116 results['dirty_image_max']['op'] = op 01117 results['dirty_image_max']['tol'] = tol 01118 01119 results['dirty_image_offsrc_max'] = {} 01120 results['dirty_image_offsrc_max']['name'] = 'Dirty image off-src rms' 01121 results['dirty_image_offsrc_max']['value'] = dirtyoffrms 01122 results['dirty_image_offsrc_max']['op'] = op 01123 results['dirty_image_offsrc_max']['tol'] = tol 01124 01125 results['clean_image_max'] = {} 01126 results['clean_image_max']['name'] = 'Clean image max' 01127 results['clean_image_max']['value'] = cleanmax 01128 results['clean_image_max']['op'] = op 01129 results['clean_image_max']['tol'] = tol 01130 01131 results['clean_image_offsrc_max'] = {} 01132 results['clean_image_offsrc_max']['name'] = 'Clean image off-src rms' 01133 results['clean_image_offsrc_max']['value'] = cleanoffrms 01134 results['clean_image_offsrc_max']['op'] = op 01135 results['clean_image_offsrc_max']['tol'] = tol 01136 01137 results['taper_image_max'] = {} 01138 results['taper_image_max']['name'] = 'Tapered image max' 01139 results['taper_image_max']['value'] = tapermax 01140 results['taper_image_max']['op'] = op 01141 results['taper_image_max']['tol'] = tol 01142 01143 results['taper_image_offsrc_max'] = {} 01144 results['taper_image_offsrc_max']['name'] = 'Tapered image off-src rms' 01145 results['taper_image_offsrc_max']['value'] = taperoffrms 01146 results['taper_image_offsrc_max']['op'] = op 01147 results['taper_image_offsrc_max']['tol'] = tol 01148 01149 results['clean_momentzero_max'] = {} 01150 results['clean_momentzero_max']['name'] = 'Moment 0 image max' 01151 results['clean_momentzero_max']['value'] = momzero_max 01152 results['clean_momentzero_max']['op'] = op 01153 results['clean_momentzero_max']['tol'] = tol 01154 01155 results['clean_momentzero_rms'] = {} 01156 results['clean_momentzero_rms']['name'] = 'Moment 0 image rms' 01157 results['clean_momentzero_rms']['value'] = momzero_rms 01158 results['clean_momentzero_rms']['op'] = op 01159 results['clean_momentzero_rms']['tol'] = tol 01160 01161 op = 'diff' 01162 tol = 0.1 01163 results['clean_momentone_median'] = {} 01164 results['clean_momentone_median']['name'] = 'Moment 1 image median' 01165 results['clean_momentone_median']['value'] = momone_median 01166 results['clean_momentone_median']['op'] = op 01167 results['clean_momentone_median']['tol'] = tol 01168 01169 results['clean_momentone_planezero'] = {} 01170 results['clean_momentone_planezero']['name'] = 'Moment 1 plane 0' 01171 results['clean_momentone_planezero']['value'] = momone_plane0 01172 results['clean_momentone_planezero']['op'] = op 01173 results['clean_momentone_planezero']['tol'] = tol 01174 01175 results['clean_momentone_planelast'] = {} 01176 results['clean_momentone_planelast']['name'] = 'Moment 1 plane 255' 01177 results['clean_momentone_planelast']['value'] = momone_plane255 01178 results['clean_momentone_planelast']['op'] = op 01179 results['clean_momentone_planelast']['tol'] = tol 01180 01181 op = 'divf' 01182 tol = 0.08 01183 results['vis_mean_cal'] = {} 01184 results['vis_mean_cal']['name'] = 'Vis mean of cal' 01185 results['vis_mean_cal']['value'] = vismean_cal 01186 results['vis_mean_cal']['op'] = op 01187 results['vis_mean_cal']['tol'] = tol 01188 01189 results['vis_mean_3c286'] = {} 01190 results['vis_mean_3c286']['name'] = 'Vis mean of 3c286' 01191 results['vis_mean_3c286']['value'] = vismean_3c286 01192 results['vis_mean_3c286']['op'] = op 01193 results['vis_mean_3c286']['tol'] = tol 01194 01195 results['model_clean_sum'] = {} 01196 results['model_clean_sum']['name'] = 'Clean Model image sum' 01197 results['model_clean_sum']['value'] = cleanmodflux 01198 results['model_clean_sum']['op'] = op 01199 results['model_clean_sum']['tol'] = tol 01200 01201 results['model_taper_sum'] = {} 01202 results['model_taper_sum']['name'] = 'Tapered Model image sum' 01203 results['model_taper_sum']['value'] = tapermodflux 01204 results['model_taper_sum']['op'] = op 01205 results['model_taper_sum']['tol'] = tol 01206 01207 results['model_pbcor_sum'] = {} 01208 results['model_pbcor_sum']['name'] = 'PBcor Model image sum' 01209 results['model_pbcor_sum']['value'] = pbcormodflux 01210 results['model_pbcor_sum']['op'] = op 01211 results['model_pbcor_sum']['tol'] = tol 01212 01213 # Now go through and regress 01214 resultlist = ['dirty_image_max','dirty_image_offsrc_max', 01215 'clean_image_max','clean_image_offsrc_max', 01216 'taper_image_max','taper_image_offsrc_max', 01217 'clean_momentzero_max','clean_momentzero_rms','clean_momentone_median', 01218 'clean_momentone_planezero','clean_momentone_planelast', 01219 'vis_mean_cal','vis_mean_3c286', 01220 'model_clean_sum','model_taper_sum','model_pbcor_sum'] 01221 01222 for keys in resultlist: 01223 res = results[keys] 01224 if prev_results.has_key(keys): 01225 # This is a known regression 01226 prev = prev_results[keys] 01227 new_val = res['value'] 01228 prev_val = prev['value'] 01229 if res['op'] == 'divf': 01230 new_diff = (new_val - prev_val)/prev_val 01231 else: 01232 new_diff = new_val - prev_val 01233 01234 if abs(new_diff)>res['tol']: 01235 new_status = 'Failed' 01236 else: 01237 new_status = 'Passed' 01238 01239 results[keys]['prev'] = prev_val 01240 results[keys]['diff'] = new_diff 01241 results[keys]['status'] = new_status 01242 results[keys]['test'] = 'Last' 01243 elif canonical_results.has_key(keys): 01244 # Go back to canonical values 01245 prev = canonical_results[keys] 01246 new_val = res['value'] 01247 prev_val = prev['value'] 01248 if res['op'] == 'divf': 01249 new_diff = (new_val - prev_val)/prev_val 01250 else: 01251 new_diff = new_val - prev_val 01252 01253 if abs(new_diff)>res['tol']: 01254 new_status = 'Failed' 01255 else: 01256 new_status = 'Passed' 01257 01258 results[keys]['prev'] = prev_val 01259 results[keys]['diff'] = new_diff 01260 results[keys]['status'] = new_status 01261 results[keys]['test'] = 'Canon' 01262 else: 01263 # Unknown regression key 01264 results[keys]['prev'] = 0.0 01265 results[keys]['diff'] = 1.0 01266 results[keys]['status'] = 'Missed' 01267 results[keys]['test'] = 'none' 01268 01269 # Done filling results 01270 new_regression['results'] = results 01271 01272 # 01273 ########################################################################## 01274 # Now the timing results 01275 ########################################################################## 01276 # 01277 datasize_raw = sdmsize 01278 datasize_fill = fillsize # (after fill) 01279 datasize_ms = fillsize 01280 01281 datasize_split = splitsize # (after split) 01282 datasize_ms = splitsize 01283 01284 datasize_clear = clearsize # (after clearcal) 01285 datasize_ms = clearsize 01286 01287 datasize_final = msallsize # (after final split) 01288 01289 datasize_fits = fitssize # (after export) 01290 01291 new_regression['datasize'] = {} 01292 new_regression['datasize']['raw'] = datasize_raw 01293 new_regression['datasize']['filled'] = datasize_fill 01294 new_regression['datasize']['ms'] = datasize_ms 01295 new_regression['datasize']['split'] = datasize_split 01296 new_regression['datasize']['clear'] = datasize_clear 01297 new_regression['datasize']['final'] = datasize_final 01298 new_regression['datasize']['fits'] = datasize_fits 01299 01300 # Save timing to regression dictionary 01301 timing = {} 01302 01303 total = {} 01304 total['wall'] = (endTime - startTime) 01305 total['cpu'] = (endProc - startProc) 01306 total['rate_raw'] = (datasize_raw/(endTime - startTime)) 01307 total['rate_ms'] = (datasize_ms/(endTime - startTime)) 01308 01309 nstages = stagetime.__len__() 01310 timing['total'] = total 01311 timing['nstages'] = nstages 01312 timing['stagename'] = stagename 01313 timing['stagetime'] = stagetime 01314 01315 new_regression['timing'] = timing 01316 01317 # 01318 ########################################################################## 01319 # Save regression results as dictionary using Pickle 01320 # 01321 pickfile = 'out.'+prefix + '.regression.'+datestring+'.pickle' 01322 f = open(pickfile,'w') 01323 p = pickle.Pickler(f) 01324 p.dump(new_regression) 01325 f.close() 01326 01327 print "" 01328 print "Regression result dictionary saved in "+pickfile 01329 print "" 01330 print "Use Pickle to retrieve these" 01331 print "" 01332 01333 # e.g. 01334 # f = open(pickfile) 01335 # u = pickle.Unpickler(f) 01336 # clnmodel = u.load() 01337 # polmodel = u.load() 01338 # f.close() 01339 01340 # 01341 ########################################################################## 01342 # Printing out results 01343 ########################################################################## 01344 # 01345 print '' 01346 print >>logfile,'' 01347 01348 # Print out comparison: 01349 res = {} 01350 resultlist = ['dirty_image_max','dirty_image_offsrc_max', 01351 'clean_image_max','clean_image_offsrc_max', 01352 'taper_image_max','taper_image_offsrc_max', 01353 'clean_momentzero_max','clean_momentzero_rms','clean_momentone_median', 01354 'clean_momentone_planezero','clean_momentone_planelast', 01355 'vis_mean_cal','vis_mean_3c286', 01356 'model_clean_sum','model_taper_sum','model_pbcor_sum'] 01357 01358 # First versus canonical values 01359 print >>logfile,'---' 01360 print >>logfile,'Regression versus previous values:' 01361 print >>logfile,'---' 01362 print '---' 01363 print 'Regression versus previous values:' 01364 print '---' 01365 01366 if regression['exist']: 01367 print >>logfile," Regression results filled from "+regressfile 01368 print >>logfile," Regression from version "+regression['version']+" on "+regression['date'] 01369 print >>logfile," Regression platform "+regression['host'] 01370 01371 print " Regression results filled from "+regressfile 01372 print " Regression from version "+regression['version']+" on "+regression['date'] 01373 print " Regression platform "+regression['host'] 01374 if regression.has_key('aipspath'): 01375 print >>logfile," Regression casapath "+regression['aipspath'] 01376 print " Regression casapath "+regression['aipspath'] 01377 01378 else: 01379 print >>logfile," No previous regression file" 01380 01381 print "" 01382 print >>logfile,"" 01383 01384 final_status = 'Passed' 01385 for keys in resultlist: 01386 res = results[keys] 01387 print '--%30s : %12.6f was %12.6f %4s %12.6f (%6s) %s ' % ( res['name'], res['value'], res['prev'], res['op'], res['diff'], res['status'], res['test'] ) 01388 print >>logfile,'--%30s : %12.6f was %12.6f %4s %12.6f (%6s) %s ' % ( res['name'], res['value'], res['prev'], res['op'], res['diff'], res['status'], res['test'] ) 01389 if res['status']=='Failed': 01390 final_status = 'Failed' 01391 01392 if (final_status == 'Passed'): 01393 regstate=True 01394 print >>logfile,'---' 01395 print >>logfile,'Passed Regression test for LeoRing' 01396 print >>logfile,'---' 01397 print 'Passed Regression test for LeoRing' 01398 else: 01399 regstate=False 01400 print >>logfile,'----FAILED Regression test for LeoRing' 01401 print '----FAILED Regression test for LeoRing' 01402 01403 print '' 01404 print 'Total wall clock time was: '+str(endTime - startTime) 01405 print 'Total CPU time was: '+str(endProc - startProc) 01406 print >>logfile,'' 01407 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 01408 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 01409 01410 print '' 01411 print 'Canonical wall clock time was: '+str(2354.63)+' on 2010-01-29 (rishi)' 01412 print 'Canonical CPU time was: '+str(2354.63)+' on 2010-01-29 (rishi)' 01413 print >>logfile,'' 01414 print >>logfile,'Canonical wall clock time was: '+str(2354.63)+' on 2010-01-29 (rishi)' 01415 print >>logfile,'Canonical CPU time was: '+str(2354.63)+' on 2010-01-29 (rishi)' 01416 01417 print '' 01418 print 'SDM datasize (MB) : '+str(datasize_raw) 01419 print 'Filled MS datasize (MB) : '+str(datasize_fill) 01420 print >>logfile,'' 01421 print >>logfile,'SDM datasize (MB) : '+str(datasize_raw) 01422 print >>logfile,'Filled MS datasize (MB) : '+str(datasize_fill) 01423 01424 print 'Split MS datasize (MB) : '+str(datasize_split) 01425 print >>logfile,'Split MS datasize (MB) : '+str(datasize_split) 01426 01427 print 'Clearcal MS datasize (MB) : '+str(datasize_clear) 01428 print >>logfile,'Clearcal MS datasize (MB) : '+str(datasize_clear) 01429 01430 print 'Final MS datasize (MB) : '+str(datasize_final) 01431 print >>logfile,'Clearcal MS datasize (MB) : '+str(datasize_final) 01432 01433 print 'Export UVFITS datasize (MB) : '+str(datasize_fits) 01434 print >>logfile,'Export UVFITS datasize (MB) : '+str(datasize_fits) 01435 01436 print 'Reported sizes are in MB from "du -ms" (1024x1024 bytes)' 01437 print >>logfile,'Reported sizes are in MB from "du -ms" (1024x1024 bytes)' 01438 01439 print '' 01440 print >>logfile,'' 01441 01442 print 'Net SDMtoMS filling I/O rate MB/s : '+str((datasize_raw+datasize_fill)/fillTime) 01443 print >>logfile,'Net SDMtoMS filling I/O rate MB/s : '+str((datasize_raw+datasize_fill)/fillTime) 01444 01445 print 'Net Splitting I/O rate MB/s was : '+str((datasize_fill+datasize_split)/splitTime) 01446 print >>logfile,'Net Splitting I/O rate MB/s was : '+str((datasize_fill+datasize_split)/splitTime) 01447 01448 print 'Net Clearcal I/O rate MB/s was : '+str((datasize_ms)/clearTime) 01449 print >>logfile,'Net Clearcal I/O rate MB/s was : '+str((datasize_ms)/clearTime) 01450 01451 print 'Net MStoUVFITS export I/O rate MB/s : '+str((datasize_final+datasize_fits)/exportTime) 01452 print >>logfile,'Net MStoUVFITS export I/O rate MB/s : '+str((datasize_final+datasize_fits)/exportTime) 01453 print '' 01454 print >>logfile,'' 01455 print 'Net SDM processing rate MB/s : '+str(datasize_raw/(endTime - startTime)) 01456 print >>logfile,'Net SDM processing rate MB/s : '+str(datasize_raw/(endTime - startTime)) 01457 print 'Net MS processing rate MB/s : '+str(datasize_ms/(endTime - startTime)) 01458 print >>logfile,'Net MS processing rate MB/s : '+str(datasize_ms/(endTime - startTime)) 01459 01460 print '* Breakdown: *' 01461 print >>logfile,'* Breakdown: *' 01462 01463 for i in range(nstages): 01464 print '* %40s * time was: %10.3f ' % (stagename[i],stagetime[i]) 01465 print >>logfile,'* %40s * time was: %10.3f ' % (stagename[i],stagetime[i]) 01466 01467 logfile.close() 01468 01469 print "Done with regression for "+mydataset 01470 # End of script 01471 # 01472 #================================================================================ 01473 #MeasurementSet Name: /home/rishi/smyers/LeoRing/leo2pt.55183.452640752315.ms 01474 #MS Version 2 01475 #================================================================================ 01476 # Observer: leo2pt Project: T.B.D. 01477 #Observation: EVLA 01478 #Data records: 745470 Total integration time = 11294 seconds 01479 # Observed from 18-Dec-2009/10:51:51.5 to 18-Dec-2009/14:00:05.5 (UTC) 01480 # 01481 # ObservationID = 0 ArrayID = 0 01482 # Date Timerange (UTC) Scan FldId FieldName nVis Int(s) SpwIds 01483 # 18-Dec-2009/10:51:51.5 - 10:53:47.5 1 0 J1042+1203 7722 1 [0] 01484 # 10:53:48.5 - 10:56:53.5 2 1 J1042+1203 12276 1 [0] 01485 # 10:56:54.5 - 11:06:59.5 3 2 Leo-1 39996 1 [0] 01486 # 11:07:00.5 - 11:17:05.5 4 3 Leo-2 39996 1 [0] 01487 # 11:17:06.5 - 11:22:12.5 5 2 Leo-1 20262 1 [0] 01488 # 11:22:13.5 - 11:27:18.5 6 3 Leo-2 20196 1 [0] 01489 # 11:27:19.5 - 11:30:26.5 7 1 J1042+1203 12408 1 [0] 01490 # 11:30:27.5 - 11:40:32.5 8 2 Leo-1 39996 1 [0] 01491 # 11:40:33.5 - 11:50:38.5 9 3 Leo-2 39996 1 [0] 01492 # 11:50:39.5 - 11:55:44.5 10 2 Leo-1 20196 1 [0] 01493 # 11:55:45.5 - 12:00:51.5 11 3 Leo-2 20262 1 [0] 01494 # 12:00:52.5 - 12:03:58.5 12 1 J1042+1203 12342 1 [0] 01495 # 12:03:59.5 - 12:14:05.5 13 2 Leo-1 40062 1 [0] 01496 # 12:14:06.5 - 12:24:10.5 14 3 Leo-2 39930 1 [0] 01497 # 12:24:11.5 - 12:29:17.5 15 2 Leo-1 20262 1 [0] 01498 # 12:29:18.5 - 12:34:23.5 16 3 Leo-2 20196 1 [0] 01499 # 12:34:24.5 - 12:37:31.5 17 1 J1042+1203 12408 1 [0] 01500 # 12:37:32.5 - 12:47:37.5 18 2 Leo-1 39996 1 [0] 01501 # 12:47:38.5 - 12:57:43.5 19 3 Leo-2 39996 1 [0] 01502 # 12:57:44.5 - 13:02:50.5 20 2 Leo-1 20262 1 [0] 01503 # 13:02:51.5 - 13:07:56.5 21 3 Leo-2 20196 1 [0] 01504 # 13:07:57.5 - 13:11:04.5 22 1 J1042+1203 12408 1 [0] 01505 # 13:11:05.5 - 13:21:10.5 23 2 Leo-1 39996 1 [0] 01506 # 13:21:11.5 - 13:31:16.5 24 3 Leo-2 39996 1 [0] 01507 # 13:31:17.5 - 13:36:23.5 25 2 Leo-1 20262 1 [0] 01508 # 13:36:24.5 - 13:41:29.5 26 3 Leo-2 20196 1 [0] 01509 # 13:41:30.5 - 13:44:37.5 27 1 J1042+1203 12408 1 [0] 01510 # 13:44:38.5 - 14:00:05.5 28 4 J1331+3030 61248 1 [0] 01511 # (nVis = Total number of time/baseline visibilities per scan) 01512 #Fields: 5 01513 # ID Code Name RA Decl Epoch SrcId nVis 01514 # 0 NONE J1042+1203 10:42:44.6052 +12.03.31.2641 J2000 0 7722 01515 # 1 D J1042+1203 10:42:44.6052 +12.03.31.2641 J2000 1 74250 01516 # 2 NONE Leo-1 10:47:22.0000 +12.16.38.0000 J2000 2 301290 01517 # 3 NONE Leo-2 10:46:45.0000 +11.50.38.0000 J2000 3 300960 01518 # 4 K J1331+3030 13:31:08.2880 +30.30.32.9589 J2000 4 61248 01519 # (nVis = Total number of time/baseline visibilities per field) 01520 #Spectral Windows: (1 unique spectral windows and 1 unique polarization setups) 01521 # SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 01522 # 0 256 TOPO 1415.3756 7.8125 2000 1415.3756 RR LL 01523 #Sources: 6 01524 # ID Name SpwId RestFreq(MHz) SysVel(km/s) 01525 # 0 J1042+1203 0 - - 01526 # 1 J1042+1203 0 - - 01527 # 2 Leo-1 0 - - 01528 # 3 Leo-2 0 - - 01529 # 4 J1331+3030 0 - - 01530 # 5 J1331+3030 0 - - 01531 #Antennas: 12: 01532 # ID Name Station Diam. Long. Lat. 01533 # 0 ea02 E02 25.0 m -107.37.04.4 +33.54.01.1 01534 # 1 ea03 E09 25.0 m -107.36.45.1 +33.53.53.6 01535 # 2 ea04 W01 25.0 m -107.37.05.9 +33.54.00.5 01536 # 3 ea05 W08 25.0 m -107.37.21.6 +33.53.53.0 01537 # 4 ea08 N01 25.0 m -107.37.06.0 +33.54.01.8 01538 # 5 ea09 E06 25.0 m -107.36.55.6 +33.53.57.7 01539 # 6 ea15 W06 25.0 m -107.37.15.6 +33.53.56.4 01540 # 7 ea19 W04 25.0 m -107.37.10.8 +33.53.59.1 01541 # 8 ea24 W05 25.0 m -107.37.13.0 +33.53.57.8 01542 # 9 ea25 N02 25.0 m -107.37.06.2 +33.54.03.5 01543 # 10 ea27 E03 25.0 m -107.37.02.8 +33.54.00.5 01544 # 11 ea28 N08 25.0 m -107.37.07.5 +33.54.15.8
1.8.0
