|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################################################## 00002 # # 00003 # Test Name: # 00004 # Use Case/Regression/Benchmark Script for NGC 5921 # 00005 # Rationale for Inclusion: # 00006 # This script is used as a running example in the CASA cookbook # 00007 # (circa 2007). It is included in the regression suite to ensure # 00008 # that modifications to CASA do not invalidate the cookbook examples. # 00009 # Features Tested: # 00010 # The script illustrates end-to-end processing with CASA # 00011 # as depicted in the following flow-chart. # 00012 # # 00013 # Input Data Process Output Data # 00014 # # 00015 # NGC5921.fits ----> importuvfits ------> ngc5921.ms + # 00016 # (1.4GHz, | ngc5921.ms.flagversions # 00017 # 63 sp chan, v # 00018 # D-array) listobs ------> casapy.log # 00019 # | # 00020 # v # 00021 # flagautocorr # 00022 # | # 00023 # v # 00024 # setjy # 00025 # | # 00026 # v # 00027 # bandpass ------> ngc5921.bcal # 00028 # | # 00029 # v # 00030 # gaincal ------> ngc5921.gcal # 00031 # | # 00032 # v # 00033 # listcal ------> ngc5921.listcal.out # 00034 # | # 00035 # v # 00036 # fluxscale ------> ngc5921.fluxscale # 00037 # | # 00038 # v # 00039 # applycal ------> ngc5921.ms # 00040 # | # 00041 # v # 00042 # split ------> ngc5921.cal.split.ms # 00043 # | # 00044 # v # 00045 # split ------> ngc5921.src.split.ms # 00046 # | # 00047 # v # 00048 # exportuvfits ------> ngc5921.split.uvfits # 00049 # | # 00050 # v # 00051 # uvcontsub ------> ngc5921.ms.cont + # 00052 # | ngc5921.ms.contsub # 00053 # v # 00054 # listvis ------> ngc5921.listvis.out # 00055 # | # 00056 # v # 00057 # clean ------> ngc5921.clean.image + # 00058 # | ngc5921.clean.model + # 00059 # | ngc5921.clean.residual # 00060 # v # 00061 # exportfits ------> ngc5921.clean.fits # 00062 # | # 00063 # v # 00064 # imstat ------> casapy.log # 00065 # | # 00066 # v # 00067 # immoments ------> ngc5921.moments.integrated + # 00068 # | ngc5921.moments.weighted_coord # 00069 # v # 00070 # Get various statistics # 00071 # on ms, image, moment images --> See benchmarking log # 00072 # Existence tests for UVFITS # 00073 # and image FITS # 00074 # Success/failure criteria: # 00075 # Listcal and Listvis output files match output from past runs PASS/FAIL # 00076 # Existence of exported UVFITS and FITS image files PASS/FAIL # 00077 # UVFITS max value test PASS/FAIL # 00078 # FITS image max and rms test PASS/FAIL # 00079 # Cal max ampl, src max ampl, image max and rms tests PASS/FAIL # 00080 # Benchmark results compared to past runs by external framework # 00081 # # 00082 ############################################################################## 00083 # # 00084 # Use Case/Regression/Benchmark Script for NGC 5921 # 00085 # # 00086 # Converted by RRusk 2007-10-31 from ngc5921_usecase.py # 00087 # Updated RRusk 2007-11-01 import regression_utility.py # 00088 # Tests loading of uvfits file # 00089 # Updated SMyers 2007-11-08 better import/export testing # 00090 # different pass/fail bookkeeping # 00091 # Updated RRusk 2007-11-08 Added some test template info # 00092 # Updated RRusk 2007-11-08 More teplate info # 00093 # Updated JCrossley 2008-09-17 Added listvis and listcal tests # 00094 # Converted by RReid 2010-01-31 from ngc5921_regression.py # 00095 # Updated RReid 2010-05-25 Made listvis use uvcontsub2 output # 00096 ############################################################################## 00097 00098 import os 00099 import time 00100 import regression_utility as tstutl 00101 import listing 00102 import datetime 00103 00104 # Enable benchmarking? 00105 benchmarking = True 00106 00107 # Run exportuvfits asynchronously (twice)? 00108 export_asynchronously = True 00109 00110 checklistvis=True 00111 # 00112 # Set up some useful variables 00113 # 00114 # Get path to CASA home directory by stipping name from '$CASAPATH' 00115 pathname=os.environ.get('CASAPATH').split()[0] 00116 00117 # This is where the NGC5921 UVFITS data will be 00118 fitsdata=pathname+'/data/regression/ngc5921/ngc5921.fits' 00119 00120 # The testdir where all output files will be kept 00121 testdir='ngc5921_regression' 00122 00123 # The prefix to use for all output files. 00124 prefix=testdir+"/"+'ngc5921' 00125 00126 # Make new test directory 00127 # (WARNING! Removes old test directory of the same name if one exists) 00128 tstutl.maketestdir(testdir) 00129 00130 datestring = datetime.datetime.isoformat(datetime.datetime.today()) 00131 outfile = 'ngc5921.' + datestring + '.log' 00132 n5921reglog = open(outfile, 'w') 00133 00134 def reportresults(redi): 00135 """ 00136 Helper function to pretty print the results in redi and report whether 00137 they all passed. redi is a dict with format 00138 {test_name1: (pass/fail, 00139 (Optional) note, 00140 (Even more optional) quantitative difference), 00141 test_name2: (pass/fail, 00142 (Optional) note, 00143 (Even more optional) quantitative difference)} 00144 """ 00145 passfail = {True: '* Passed', False: '* FAILED'} 00146 normalsevere = {True: 'NORMAL', False: 'SEVERE'} 00147 ok = True 00148 00149 for t in redi: 00150 tup = redi[t] 00151 msg = passfail[tup[0]] + ' ' + t + ' test' 00152 print >>n5921reglog, msg 00153 tstutl.note(msg, normalsevere[tup[0]]) 00154 if len(tup) > 1: 00155 print >>n5921reglog, tup[1] 00156 #tstutl.note("\"tup[0]\": \"%s\"" % tup[0], "WARN") 00157 if not tup[0]: 00158 ok = False 00159 if len(tup) > 2: 00160 tstutl.note(' ' + t + " difference: " + str(tup[2]), 00161 normalsevere[tup[0]]) 00162 return ok 00163 00164 def listfailures(redi): 00165 """ 00166 Helper function to summarize any failures in redi at the end. 00167 redi has the same format as in reportresults. 00168 """ 00169 for t in redi: 00170 tup = redi[t] 00171 if not tup[0]: 00172 msg = t + " FAILED" 00173 if len(tup) > 2: 00174 msg += ":\n " + str(tup[1]) + "\n difference: " + str(tup[2]) 00175 tstutl.note(msg, "SEVERE") 00176 00177 # Start benchmarking 00178 if benchmarking: 00179 startTime = time.time() 00180 startProc = time.clock() 00181 00182 passedall = True # So far! 00183 00184 # 00185 #===================================================================== 00186 # 00187 # Import the data from FITS to MS 00188 # 00189 print '--Import--' 00190 00191 # Safest to start from task defaults 00192 default('importuvfits') 00193 00194 # Set up the MS filename and save as new global variable 00195 msfile = prefix + '.ms' 00196 00197 # Use task importuvfits 00198 fitsfile = fitsdata 00199 vis = msfile 00200 antnamescheme = 'new' 00201 importuvfits() 00202 00203 # Note that there will be a ngc5921.ms.flagversions 00204 # containing the initial flags as backup for the main ms flags. 00205 00206 # Record import time 00207 if benchmarking: 00208 importtime = time.time() 00209 00210 # 00211 #===================================================================== 00212 # 00213 # List a summary of the MS 00214 # 00215 print '--Listobs--' 00216 00217 # Don't default this one. Make use of the previous setting of 00218 # vis. Remember, the variables are GLOBAL! 00219 00220 # You may wish to see more detailed information, like the scans. 00221 # In this case use the verbose = True option 00222 verbose = True 00223 00224 listobs() 00225 00226 # You should get in your logger window and in the casapy.log file 00227 # something like: 00228 # 00229 # MeasurementSet Name: /home/sandrock2/smyers/Testing2/Sep07/ngc5921_regression/ngc5921.ms 00230 # MS Version 2 00231 # 00232 # Observer: TEST Project: 00233 # Observation: VLA 00234 # 00235 # Data records: 22653 Total integration time = 5280 seconds 00236 # Observed from 09:19:00 to 10:47:00 00237 # 00238 # ObservationID = 0 ArrayID = 0 00239 # Date Timerange Scan FldId FieldName SpwIds 00240 # 13-Apr-1995/09:19:00.0 - 09:24:30.0 1 0 1331+30500002_0 [0] 00241 # 09:27:30.0 - 09:29:30.0 2 1 1445+09900002_0 [0] 00242 # 09:33:00.0 - 09:48:00.0 3 2 N5921_2 [0] 00243 # 09:50:30.0 - 09:51:00.0 4 1 1445+09900002_0 [0] 00244 # 10:22:00.0 - 10:23:00.0 5 1 1445+09900002_0 [0] 00245 # 10:26:00.0 - 10:43:00.0 6 2 N5921_2 [0] 00246 # 10:45:30.0 - 10:47:00.0 7 1 1445+09900002_0 [0] 00247 # 00248 # Fields: 3 00249 # ID Code Name Right Ascension Declination Epoch 00250 # 0 C 1331+30500002_013:31:08.29 +30.30.32.96 J2000 00251 # 1 A 1445+09900002_014:45:16.47 +09.58.36.07 J2000 00252 # 2 N5921_2 15:22:00.00 +05.04.00.00 J2000 00253 # 00254 # Spectral Windows: (1 unique spectral windows and 1 unique polarization setups) 00255 # SpwID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 00256 # 0 63 LSRK 1412.68608 24.4140625 1550.19688 1413.44902 RR LL 00257 # 00258 # Feeds: 28: printing first row only 00259 # Antenna Spectral Window # Receptors Polarizations 00260 # 1 -1 2 [ R, L] 00261 # 00262 # Antennas: 27: 00263 # ID Name Station Diam. Long. Lat. 00264 # 0 1 VLA:N7 25.0 m -107.37.07.2 +33.54.12.9 00265 # 1 2 VLA:W1 25.0 m -107.37.05.9 +33.54.00.5 00266 # 2 3 VLA:W2 25.0 m -107.37.07.4 +33.54.00.9 00267 # 3 4 VLA:E1 25.0 m -107.37.05.7 +33.53.59.2 00268 # 4 5 VLA:E3 25.0 m -107.37.02.8 +33.54.00.5 00269 # 5 6 VLA:E9 25.0 m -107.36.45.1 +33.53.53.6 00270 # 6 7 VLA:E6 25.0 m -107.36.55.6 +33.53.57.7 00271 # 7 8 VLA:W8 25.0 m -107.37.21.6 +33.53.53.0 00272 # 8 9 VLA:N5 25.0 m -107.37.06.7 +33.54.08.0 00273 # 9 10 VLA:W3 25.0 m -107.37.08.9 +33.54.00.1 00274 # 10 11 VLA:N4 25.0 m -107.37.06.5 +33.54.06.1 00275 # 11 12 VLA:W5 25.0 m -107.37.13.0 +33.53.57.8 00276 # 12 13 VLA:N3 25.0 m -107.37.06.3 +33.54.04.8 00277 # 13 14 VLA:N1 25.0 m -107.37.06.0 +33.54.01.8 00278 # 14 15 VLA:N2 25.0 m -107.37.06.2 +33.54.03.5 00279 # 15 16 VLA:E7 25.0 m -107.36.52.4 +33.53.56.5 00280 # 16 17 VLA:E8 25.0 m -107.36.48.9 +33.53.55.1 00281 # 17 18 VLA:W4 25.0 m -107.37.10.8 +33.53.59.1 00282 # 18 19 VLA:E5 25.0 m -107.36.58.4 +33.53.58.8 00283 # 19 20 VLA:W9 25.0 m -107.37.25.1 +33.53.51.0 00284 # 20 21 VLA:W6 25.0 m -107.37.15.6 +33.53.56.4 00285 # 21 22 VLA:E4 25.0 m -107.37.00.8 +33.53.59.7 00286 # 23 24 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1 00287 # 24 25 VLA:N6 25.0 m -107.37.06.9 +33.54.10.3 00288 # 25 26 VLA:N9 25.0 m -107.37.07.8 +33.54.19.0 00289 # 26 27 VLA:N8 25.0 m -107.37.07.5 +33.54.15.8 00290 # 27 28 VLA:W7 25.0 m -107.37.18.4 +33.53.54.8 00291 # 00292 # Tables: 00293 # MAIN 22653 rows 00294 # ANTENNA 28 rows 00295 # DATA_DESCRIPTION 1 row 00296 # DOPPLER <absent> 00297 # FEED 28 rows 00298 # FIELD 3 rows 00299 # FLAG_CMD <empty> 00300 # FREQ_OFFSET <absent> 00301 # HISTORY 273 rows 00302 # OBSERVATION 1 row 00303 # POINTING 168 rows 00304 # POLARIZATION 1 row 00305 # PROCESSOR <empty> 00306 # SOURCE 3 rows 00307 # SPECTRAL_WINDOW 1 row 00308 # STATE <empty> 00309 # SYSCAL <absent> 00310 # WEATHER <absent> 00311 # 00312 00313 # Record listing completion time 00314 if benchmarking: 00315 listtime = time.time() 00316 00317 # 00318 #===================================================================== 00319 # 00320 # Get rid of the autocorrelations from the MS 00321 # 00322 #print '--Flagautocorr--' 00323 print '--Flag auto-correlations--' 00324 00325 # Don't default this one either, there is only one parameter (vis) 00326 default(tflagdata) 00327 vis = msfile 00328 autocorr = True 00329 flagbackup = False 00330 tflagdata() 00331 00332 # Record flagging completion time 00333 if benchmarking: 00334 flagtime = time.time() 00335 00336 # 00337 #===================================================================== 00338 # 00339 # Set the fluxes of the primary calibrator(s) 00340 # 00341 print '--Setjy--' 00342 default('setjy') 00343 00344 vis = msfile 00345 00346 # 00347 # 1331+305 = 3C286 is our primary calibrator 00348 # Use the wildcard on the end of the source name 00349 # since the field names in the MS have inherited the 00350 # AIPS qualifiers 00351 field = '1331+305*' 00352 00353 # This is 1.4GHz D-config and 1331+305 is sufficiently unresolved 00354 # that we do not need a model image. For higher frequencies 00355 # (particularly in A and B config) you would want to use one. 00356 modimage = '' 00357 00358 # Setjy knows about this source so we don't need anything more 00359 standard='Perley-Taylor 99' # enforce older standard 00360 00361 scalebychan=False 00362 00363 usescratch=False 00364 00365 setjy() 00366 00367 # 00368 # You should see something like this in the logger and casapy.log file: 00369 # 00370 # 1331+30500002_0 spwid= 0 [I=14.76, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00371 # 00372 # So its using 14.76Jy as the flux of 1331+305 in the single Spectral Window 00373 # in this MS. 00374 00375 # Record setjy completion time 00376 if benchmarking: 00377 setjytime = time.time() 00378 00379 # 00380 #===================================================================== 00381 # 00382 # Bandpass calibration 00383 # 00384 print '--Bandpass--' 00385 default('bandpass') 00386 00387 # We can first do the bandpass on the single 5min scan on 1331+305 00388 # At 1.4GHz phase stablility should be sufficient to do this without 00389 # a first (rough) gain calibration. This will give us the relative 00390 # antenna gain as a function of frequency. 00391 00392 vis = msfile 00393 00394 # set the name for the output bandpass caltable 00395 btable = prefix + '.bcal' 00396 caltable = btable 00397 00398 # No gain tables yet 00399 gaintable = '' 00400 gainfield = '' 00401 interp = '' 00402 00403 # Use flux calibrator 1331+305 = 3C286 (FIELD_ID 0) as bandpass calibrator 00404 field = '0' 00405 # all channels 00406 spw = '' 00407 # No other selection 00408 selectdata = False 00409 00410 # In this band we do not need a-priori corrections for 00411 # antenna gain-elevation curve or atmospheric opacity 00412 # (at 8GHz and above you would want these) 00413 gaincurve = False 00414 opacity = 0.0 00415 00416 # Choose bandpass solution type 00417 # Pick standard time-binned B (rather than BPOLY) 00418 bandtype = 'B' 00419 00420 # set solution interval arbitrarily long (get single bpass) 00421 solint = 'inf' 00422 combine='scan' 00423 00424 # reference antenna Name 15 (15=VLA:N2) (Id 14) 00425 refant = 'VA15' 00426 00427 bandpass() 00428 00429 # You can use plotcal to examine the solutions 00430 #default('plotcal') 00431 #caltable = btable 00432 #yaxis = 'amp' 00433 #field = '0' 00434 #overplot = True 00435 #plotcal() 00436 # 00437 #yaxis = 'phase' 00438 #plotcal() 00439 # 00440 # Note the rolloff in the start and end channels. Looks like 00441 # channels 6-56 (out of 0-62) are the best 00442 00443 # bandpass calibration completion time 00444 if benchmarking: 00445 bptime = time.time() 00446 00447 # 00448 #===================================================================== 00449 # 00450 # Gain calibration 00451 # 00452 print '--Gaincal--' 00453 default('gaincal') 00454 00455 # Armed with the bandpass, we now solve for the 00456 # time-dependent antenna gains 00457 00458 vis = msfile 00459 00460 # set the name for the output gain caltable 00461 gtable = prefix + '.gcal' 00462 caltable = gtable 00463 00464 # Use our previously determined bandpass 00465 # Note this will automatically be applied to all sources 00466 # not just the one used to determine the bandpass 00467 gaintable = btable 00468 gainfield = '0' 00469 00470 # Use nearest (there is only one bandpass entry) 00471 interp = 'nearest' 00472 00473 # Gain calibrators are 1331+305 and 1445+099 (FIELD_ID 0 and 1) 00474 field = '0,1' 00475 00476 # We have only a single spectral window (SPW 0) 00477 # Choose 51 channels 6-56 out of the 63 00478 # to avoid end effects. 00479 # Channel selection is done inside spw 00480 spw = '0:6~56' 00481 00482 # No other selection 00483 selectdata = False 00484 00485 # In this band we do not need a-priori corrections for 00486 # antenna gain-elevation curve or atmospheric opacity 00487 # (at 8GHz and above you would want these) 00488 gaincurve = False 00489 opacity = 0.0 00490 00491 # scan-based G solutions for both amplitude and phase 00492 gaintype = 'G' 00493 solint = 'inf' 00494 combine='' 00495 calmode = 'ap' 00496 00497 # minimum SNR allowed 00498 minsnr = 1.0 00499 00500 # reference antenna 15 (15=VLA:N2) 00501 refant = '15' 00502 00503 gaincal() 00504 00505 # You can use plotcal to examine the gain solutions 00506 #default('plotcal') 00507 #caltable = gtable 00508 #yaxis = 'amp' 00509 #field = '0,1' 00510 #overplot = True 00511 #plotcal() 00512 # 00513 #yaxis = 'phase' 00514 #plotcal() 00515 # 00516 # The amp and phase coherence looks good 00517 00518 # gaincal calibration completion time 00519 if benchmarking: 00520 gaintime = time.time() 00521 00522 # 00523 #===================================================================== 00524 # List calibration solutions 00525 # 00526 print '--Listcal--' 00527 listcalOut = prefix + '.listcal.out' 00528 00529 default('listcal') 00530 vis = msfile 00531 caltable = gtable 00532 listfile = listcalOut 00533 listcal() 00534 00535 # Record calibration listing time 00536 if benchmarking: 00537 listcaltime = time.time() 00538 00539 # Before testing listcal output, remove first line of file 00540 # (First line contains hard-coded path to input files) 00541 os.system('mv ' + listcalOut + ' ' + listcalOut + '.tmp') 00542 os.system('tail -n +2 ' + listcalOut + '.tmp > ' + listcalOut) 00543 os.system('rm -f ' + listcalOut + '.tmp') 00544 00545 # Test the listcal output 00546 print "Comparing listcal output with standard..." 00547 standardOut = pathname+'/data/regression/ngc5921/listcal.default.out' 00548 listcalresults = {} 00549 try: 00550 print " 1. Checking that metadata agree..." 00551 listcalresults['listcal metadata'] = (listing.diffMetadata(listcalOut, 00552 standardOut, 00553 prefix=prefix + ".listcal"),) 00554 # Test data (floats) 00555 print " 2. Checking that data agree to within allowed imprecision..." 00556 precision = '0.003' 00557 print " Allowed visibility imprecision is " + precision 00558 listcalresults['listcal data'] = (listing.diffAmpPhsFloat(listcalOut, 00559 standardOut, 00560 prefix = prefix+".listcal", 00561 precision = precision),) 00562 passedall = reportresults(listcalresults) and passedall 00563 except Exception, e: 00564 print "Error", e, "checking listcal." 00565 raise e 00566 00567 # 00568 #===================================================================== 00569 # 00570 # Bootstrap flux scale 00571 # 00572 print '--Fluxscale--' 00573 default('fluxscale') 00574 00575 vis = msfile 00576 00577 # set the name for the output rescaled caltable 00578 ftable = prefix + '.fluxscale' 00579 fluxtable = ftable 00580 00581 # point to our first gain cal table 00582 caltable = gtable 00583 00584 # we will be using 1331+305 (the source we did setjy on) as 00585 # our flux standard reference - note its extended name as in 00586 # the FIELD table summary above (it has a VLA seq number appended) 00587 reference = '1331*' 00588 00589 # we want to transfer the flux to our other gain cal source 1445+099 00590 transfer = '1445*' 00591 00592 # testing new option 00593 incremental=True 00594 00595 fluxscale() 00596 00597 # In the logger you should see something like: 00598 # Flux density for 1445+09900002_0 in SpW=0 is: 00599 # 2.48576 +/- 0.00123122 (SNR = 2018.94, nAnt= 27) 00600 00601 # If you run plotcal() on the caltable = 'ngc5921.fluxscale' 00602 # you will see now it has brought the amplitudes in line between 00603 # the first scan on 1331+305 and the others on 1445+099 00604 00605 # Record fluxscale completion time 00606 if benchmarking: 00607 fstime = time.time() 00608 00609 # 00610 #===================================================================== 00611 # 00612 # Apply our calibration solutions to the data 00613 # (This will put calibrated data into the CORRECTED_DATA column) 00614 # 00615 print '--ApplyCal--' 00616 default('applycal') 00617 00618 vis = msfile 00619 00620 # We want to correct the calibrators using themselves 00621 # and transfer from 1445+099 to itself and the target N5921 00622 00623 # Start with the fluxscale/gain and bandpass tables 00624 #gaintable = [ftable,btable] 00625 #for incremental fluxscale table case 00626 gaintable = [ftable,gtable,btable] 00627 00628 # use nearest (on sky) for gain transfer 00629 # use explicit selection on the bandpass table (all of it, in this case) 00630 #gainfield = ['nearest','0'] 00631 gainfield = ['nearest','nearest','0'] 00632 00633 # interpolation using linear for gain, nearest for bandpass 00634 #interp = ['linear','nearest'] 00635 interp = ['linear', 'linear','nearest'] 00636 00637 # only one spw, do not need mapping 00638 spwmap = [] 00639 00640 # all channels 00641 spw = '' 00642 selectdata = False 00643 00644 # as before 00645 gaincurve = False 00646 opacity = 0.0 00647 00648 # select all fields 00649 field = '' 00650 00651 applycal() 00652 00653 # Record applycal completion time 00654 if benchmarking: 00655 correcttime = time.time() 00656 00657 # 00658 #===================================================================== 00659 # 00660 # Split the gain calibrater data, then the target 00661 # 00662 print '--Split 1445+099 Data--' 00663 default('split') 00664 00665 vis = msfile 00666 00667 # We first want to write out the corrected data for the calibrator 00668 00669 # Make an output vis file 00670 calsplitms = prefix + '.cal.split.ms' 00671 outputvis = calsplitms 00672 00673 # Select the 1445+099 field, all chans 00674 field = '1445*' 00675 spw = '' 00676 00677 # pick off the CORRECTED_DATA column 00678 datacolumn = 'corrected' 00679 00680 split() 00681 00682 # Record split cal completion time 00683 if benchmarking: 00684 splitcaltime = time.time() 00685 00686 # 00687 # Now split NGC5921 data (before continuum subtraction) 00688 # 00689 print '--Split NGC5921 Data--' 00690 00691 splitms = prefix + '.src.split.ms' 00692 outputvis = splitms 00693 00694 # Pick off N5921 00695 field = 'N5921*' 00696 00697 split() 00698 00699 # Record split src data completion time 00700 if benchmarking: 00701 splitsrctime = time.time() 00702 00703 # 00704 #===================================================================== 00705 # 00706 # Export the NGC5921 data as UVFITS 00707 # Start with the split file. 00708 # 00709 print '--Export UVFITS--' 00710 default('exportuvfits') 00711 00712 srcuvfits = prefix + '.split.uvfits' 00713 00714 vis = splitms 00715 fitsfile = srcuvfits 00716 00717 # Since this is a split dataset, the calibrated data is 00718 # in the DATA column already. 00719 datacolumn = 'data' 00720 00721 # Write as a multisource UVFITS (with SU table) 00722 # even though it will have only one field in it 00723 multisource = True 00724 00725 # Run asynchronously so as not to interfere with other tasks 00726 # (BETA: also avoids crash on next importuvfits) 00727 #async = True 00728 async = export_asynchronously 00729 00730 async_exportuvfits_id = exportuvfits() 00731 print "async_exportuvfits_id =", async_exportuvfits_id 00732 00733 # Record exportuvfits completion time 00734 # NOTE: If async=true this can't be used to time exportuvfits 00735 # but it is still needed as a start time for uvcontsub. 00736 if benchmarking: 00737 exportuvfitstime = time.time() 00738 00739 # 00740 #===================================================================== 00741 # 00742 # UV-plane continuum subtraction on the target 00743 # (this will update the CORRECTED_DATA column) 00744 # 00745 print '--UV Continuum Subtract--' 00746 default('uvcontsub') 00747 00748 vis = msfile 00749 00750 # Pick off N5921 00751 field = 'N5921*' 00752 00753 # Use channels 4-6 and 50-59 for continuum fit 00754 fitspw='0:4~6;50~59' 00755 00756 # Subtr and write out all of spw=0 00757 spw = '0' 00758 00759 # Averaging time (per integration) 00760 solint = 'int' 00761 00762 # Fit only a mean level 00763 fitorder = 0 00764 00765 # Let it split out the data automatically for us 00766 want_cont = True 00767 00768 uvcontsub() 00769 00770 # You will see it made two new MS: 00771 # ngc5921.ms.cont (Continuum estimate) 00772 # ngc5921.ms.contsub (Continuum subtracted) 00773 00774 srcsplitms = msfile + '.contsub' 00775 00776 # Record continuum subtraction time 00777 if benchmarking: 00778 contsubtime = time.time() 00779 00780 if checklistvis: 00781 #===================================================================== 00782 # List corrected data in MS 00783 # 00784 print '--Listvis--' 00785 listvisOut = prefix + '.listvis.out' 00786 00787 default('listvis') 00788 # vis = srcsplitms 00789 # datacolumn = 'data' 00790 # selectdata=True 00791 # antenna='VA03&VA04' 00792 # listfile=listvisOut 00793 print "Listing corrected data." 00794 print "Reducing output by selecting only baseline 3&4." 00795 # listvis() 00796 listvis(vis=srcsplitms,datacolumn='data',selectdata=True,antenna='VA03&VA04',listfile=listvisOut) 00797 00798 # Record visibility listing time 00799 if benchmarking: 00800 listvistime = time.time() 00801 00802 # Test the listvis output 00803 print "Comparing continuum subtracted listvis output with repository standard..." 00804 standardOut = pathname+'/data/regression/ngc5921/listvis.ant34.contsub.out' 00805 00806 # Test metadata 00807 print " Checking that metadata agree and that data agree to within allowed imprecision..." 00808 precision = '0.200' 00809 print " Allowed visibility imprecision is ", precision 00810 listvisresults = {} 00811 try: 00812 listvisresults['listvis metadata'] = (listing.diffMetadata(listvisOut, 00813 standardOut, 00814 prefix=prefix + ".listvis"),) 00815 listvisresults['listvis data'] = (listing.diffAmpPhsFloat(listvisOut, 00816 standardOut, 00817 prefix=prefix + ".listvis", 00818 precision=precision),) 00819 passedall = reportresults(listvisresults) and passedall 00820 except Exception, e: 00821 print "Error", e, "checking listvis." 00822 raise e 00823 00824 #===================================================================== 00825 # 00826 # Done with calibration 00827 # Now clean an image cube of N5921 00828 # 00829 print '--Clean--' 00830 default('clean') 00831 00832 # Pick up our split source data 00833 vis = srcsplitms 00834 00835 # Make an image root file name 00836 imname = prefix + '.clean' 00837 imagename = imname 00838 00839 # Set up the output image cube 00840 mode = 'channel' 00841 nchan = 46 00842 start = 5 00843 width = 1 00844 00845 # This is a single-source MS with one spw 00846 field = '0' 00847 spw = '' 00848 00849 # Standard gain factor 0.1 00850 gain = 0.1 00851 00852 # Set the output image size and cell size (arcsec) 00853 #imsize = [256,256] 00854 00855 # Do a simple Clark clean 00856 psfmode = 'clark' 00857 imagermode = '' 00858 # If desired, you can do a Cotton-Schwab clean 00859 # but will have only marginal improvement for this data 00860 #alg = 'csclean' 00861 # Twice as big for Cotton-Schwab (cleans inner quarter) 00862 #imsize = [512,512] 00863 00864 # Pixel size 15 arcsec for this data (1/3 of 45" beam) 00865 # VLA D-config L-band 00866 cell = [15.,15.] 00867 00868 # Fix maximum number of iterations 00869 niter = 6000 00870 00871 # Also set flux residual threshold (in mJy) 00872 threshold=8.0 00873 00874 # Set up the weighting 00875 # Use Briggs weighting (a moderate value, on the uniform side) 00876 weighting = 'briggs' 00877 #rmode = 'norm' 00878 robust = 0.5 00879 00880 # No clean mask or cleanbox for now 00881 mask = '' 00882 #cleanbox = [] 00883 00884 # But if you had a cleanbox saved in a file, e.g. "regionfile.txt" 00885 # you could use it: 00886 #cleanbox='regionfile.txt' 00887 # 00888 # and if you wanted to use interactive clean 00889 #cleanbox='interactive' 00890 00891 usescratch=False 00892 00893 clean() 00894 00895 # Should find stuff in the logger like: 00896 # 00897 # Fitted beam used in restoration: 51.5643 by 45.6021 (arcsec) 00898 # at pa 14.5411 (deg) 00899 # 00900 # It will have made the images: 00901 # ----------------------------- 00902 # ngc5921.clean.image 00903 # ngc5921.clean.model 00904 # ngc5921.clean.residual 00905 00906 clnimage = imname+'.image' 00907 00908 # Record clean completion time 00909 if benchmarking: 00910 cleantime = time.time() 00911 00912 # 00913 #===================================================================== 00914 # 00915 # Done with imaging 00916 # Now view the image cube of N5921 00917 # 00918 #print '--View image--' 00919 #viewer(clnimage,'image') 00920 00921 #===================================================================== 00922 # 00923 # Export the Final CLEAN Image as FITS 00924 # 00925 print '--Final Export CLEAN FITS--' 00926 default('exportfits') 00927 00928 clnfits = prefix + '.clean.fits' 00929 00930 imagename = clnimage 00931 fitsimage = clnfits 00932 00933 # Run asynchronously so as not to interfere with other tasks 00934 # (BETA: also avoids crash on next importfits) 00935 #async = True 00936 async = export_asynchronously 00937 00938 async_clean_id = exportfits() 00939 00940 # Record export FITS image completion time 00941 # NOTE: this is currently irrelevant because its async 00942 if benchmarking: 00943 exportimtime = time.time() 00944 00945 # 00946 #===================================================================== 00947 # 00948 # Get some image statistics 00949 # 00950 print '--Imstat--' 00951 default('imstat') 00952 00953 imagename = clnimage 00954 00955 cubestats =imstat() 00956 00957 # A summary of the cube will be seen in the logger 00958 # cubestats will contain a dictionary of the statistics 00959 00960 # 00961 #===================================================================== 00962 # 00963 # Get some image moments 00964 # 00965 print '--ImMoments--' 00966 default('immoments') 00967 00968 imagename = clnimage 00969 00970 # Do first and second moments 00971 moments = [0,1] 00972 00973 # Need to mask out noisy pixels, currently done 00974 # using hard global limits 00975 excludepix = [-100,0.009] 00976 00977 # Include all planes 00978 chans = '' 00979 00980 # Output root name 00981 momfile = prefix + '.moments' 00982 outfile = momfile 00983 00984 immoments() 00985 00986 momzeroimage = momfile + '.integrated' 00987 momoneimage = momfile + '.weighted_coord' 00988 00989 # It will have made the images: 00990 # -------------------------------------- 00991 # ngc5921.moments.integrated 00992 # ngc5921.moments.weighted_coord 00993 00994 # 00995 #===================================================================== 00996 # 00997 # Get some statistics of the moment images 00998 # 00999 print '--Imstat--' 01000 default('imstat') 01001 01002 imagename = momzeroimage 01003 momzerostats = imstat() 01004 01005 imagename = momoneimage 01006 momonestats = imstat() 01007 01008 # 01009 #===================================================================== 01010 # 01011 # Can do some image statistics if you wish 01012 # Treat this like a regression script 01013 # WARNING: currently requires toolkit 01014 # 01015 print "" 01016 print ' NGC5921 results ' 01017 print ' =============== ' 01018 01019 print '' 01020 print ' --Regression Tests--' 01021 print '' 01022 # 01023 # Use the ms tool to get max of the MSs 01024 # Eventually should be available from a task 01025 # 01026 # Pull the max cal amp value out of the MS 01027 ms.open(calsplitms) 01028 thistest_cal = max(ms.range(["amplitude"]).get('amplitude')) 01029 ms.close() 01030 oldtest_cal = 34.0338668823 01031 print ' Cal Max amplitude should be ',oldtest_cal 01032 print ' Found : Max = ',thistest_cal 01033 diff_cal = abs((oldtest_cal-thistest_cal)/oldtest_cal) 01034 print ' Difference (fractional) = ',diff_cal 01035 01036 print '' 01037 # Pull the max src amp value out of the MS 01038 ms.open(srcsplitms) 01039 thistest_src = max(ms.range(["amplitude"]).get('amplitude')) 01040 ms.close() 01041 #oldtest_src = 1.37963354588 # This was in chans 5-50 01042 oldtest_src = 46.2060050964 # now in all chans 01043 print ' Src Max amplitude should be ',oldtest_src 01044 print ' Found : Max = ',thistest_src 01045 diff_src = abs((oldtest_src-thistest_src)/oldtest_src) 01046 print ' Difference (fractional) = ',diff_src 01047 01048 # 01049 # Now use the stats produced by imstat above 01050 # 01051 print '' 01052 # Pull the max from the cubestats dictionary 01053 # created above using imstat 01054 thistest_immax=cubestats['max'][0] 01055 oldtest_immax = 0.052414759993553162 01056 print ' Clean image max should be ',oldtest_immax 01057 print ' Found : Image Max = ',thistest_immax 01058 diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax) 01059 print ' Difference (fractional) = ',diff_immax 01060 01061 print '' 01062 # Pull the rms from the cubestats dictionary 01063 thistest_imrms=cubestats['rms'][0] 01064 oldtest_imrms = 0.0020064469 01065 print ' Clean image rms should be ',oldtest_imrms 01066 print ' Found : Image rms = ',thistest_imrms 01067 diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms) 01068 print ' Difference (fractional) = ',diff_imrms 01069 01070 print '' 01071 # Pull the max from the momzerostats dictionary 01072 thistest_momzeromax=momzerostats['max'][0] 01073 oldtest_momzeromax = 1.4868 01074 print ' Moment 0 image max should be ',oldtest_momzeromax 01075 print ' Found : Moment 0 Max = ',thistest_momzeromax 01076 diff_momzeromax = abs(1.0 - thistest_momzeromax / oldtest_momzeromax) 01077 print ' Difference (fractional) = ',diff_momzeromax 01078 01079 print '' 01080 # Pull the mean from the momonestats dictionary 01081 thistest_momoneavg=momonestats['mean'][0] 01082 oldtest_momoneavg = 1486.8473 01083 print ' Moment 1 image mean should be ',oldtest_momoneavg 01084 print ' Found : Moment 1 Mean = ',thistest_momoneavg 01085 diff_momoneavg = abs((oldtest_momoneavg-thistest_momoneavg)/oldtest_momoneavg) 01086 print ' Difference (fractional) = ',diff_momoneavg 01087 01088 01089 # Record processing completion time 01090 if benchmarking: 01091 endProc = time.clock() 01092 endTime = time.time() 01093 01094 # 01095 #===================================================================== 01096 # Export tests 01097 # 01098 print '' 01099 print ' --Export Tests--' 01100 print '' 01101 # 01102 # First the UVFITS 01103 # 01104 tstutl.note("Opening UVFITS file " + srcuvfits + 01105 " to verify its existence") 01106 01107 if export_asynchronously: 01108 while True: 01109 try: 01110 if tm.retrieve(async_exportuvfits_id): 01111 break 01112 else: 01113 time.sleep(1) 01114 except Exception, e: 01115 tstutl.note("Error checking whether exportuvfits finished.", 01116 "SEVERE") 01117 tstutl.note("async_exportuvfits_id was " + str(async_exportuvfits_id), "SEVERE") 01118 raise e 01119 01120 uvfitsexists=False 01121 if os.path.isfile(srcuvfits): 01122 uvfitsexists=True 01123 01124 diff_uvfits = 1.0 01125 if uvfitsexists: 01126 tstutl.note("Opening the UVFITS file", "NORMAL") 01127 try: 01128 ms.fromfits(prefix+".fromuvfits.ms", srcuvfits) 01129 ms.summary() 01130 # Pull the max src amp value out of the MS 01131 fitstest_src = max(ms.range(["amplitude"]).get('amplitude')) 01132 ms.close() 01133 tstutl.note("Verified that a valid UVFITS file was written") 01134 print '' 01135 # Test the max in MS 01136 # oldtest_src = 46.2060050964 01137 print ' UVFITS Src Max amplitude should be ', thistest_src 01138 print ' Found : UVFITS Max = ', fitstest_src 01139 diff_uvfits = abs((fitstest_src - thistest_src) / thistest_src) 01140 print ' Difference (fractional) = ', diff_uvfits 01141 except Exception, e: 01142 tstutl.note("Failed to open UVFITS file because: "+e, "SEVERE") 01143 else: 01144 tstutl.note("Could not open the UVFITS file", "SEVERE") 01145 01146 # 01147 # Now the Clean FITS 01148 # 01149 print '' 01150 tstutl.note("Opening FITS image file " + fitsimage + 01151 " to verify its existence") 01152 01153 if export_asynchronously: 01154 while True: 01155 if tm.retrieve(async_clean_id): 01156 break 01157 else: 01158 time.sleep(1) 01159 01160 fitsimageexists=False 01161 if os.path.isfile(fitsimage): 01162 fitsimageexists=True 01163 01164 diff_fitsmax = 1.0 01165 diff_fitsrms = 1.0 01166 if fitsimageexists: 01167 tstutl.note("Opening the FITS image", "NORMAL") 01168 try: 01169 default('imstat') 01170 imagename = fitsimage 01171 fitstats = imstat() 01172 tstutl.note("Verified that a valid FITS image file was written") 01173 print '' 01174 # Pull the max from the fitstats dictionary 01175 fitstest_immax=fitstats['max'][0] 01176 print ' FITS image max should be ',thistest_immax 01177 print ' Found : FITS Image Max = ',fitstest_immax 01178 diff_fitsmax = abs((fitstest_immax-thistest_immax)/thistest_immax) 01179 print ' Difference (fractional) = ',diff_fitsmax 01180 print '' 01181 # Pull the rms from the fitstats dictionary 01182 fitstest_imrms=fitstats['rms'][0] 01183 print ' FITS image rms should be ',thistest_imrms 01184 print ' Found : FITS Image rms = ',fitstest_imrms 01185 diff_fitsrms = abs((fitstest_imrms-thistest_imrms)/thistest_imrms) 01186 print ' Difference (fractional) = ',diff_fitsrms 01187 except Exception, e: 01188 tstutl.note("Failed to open FITS image file because: "+e,"SEVERE") 01189 else: 01190 tstutl.note("Could not open the FITS image", "SEVERE") 01191 01192 # 01193 #===================================================================== 01194 # 01195 if not benchmarking: 01196 print '' 01197 print '--- Done ---' 01198 else: 01199 print >>n5921reglog,'' 01200 print >>n5921reglog,'' 01201 print >>n5921reglog,'********** Data Summary *********' 01202 print >>n5921reglog,'* Observer: TEST Project: *' 01203 print >>n5921reglog,'* Observation: VLA(28 antennas) *' 01204 print >>n5921reglog,'* Data records: 22653 Total integration time = 5280 seconds *' 01205 print >>n5921reglog,'* Observed from 09:19:00 to 10:47:00 *' 01206 print >>n5921reglog,'* Fields: 3 *' 01207 print >>n5921reglog,'* ID Name Right Ascension Declination Epoch *' 01208 print >>n5921reglog,'* 0 1331+30500002_013:31:08.29 +30.30.32.96 J2000 *' 01209 print >>n5921reglog,'* 1 1445+09900002_014:45:16.47 +09.58.36.07 J2000 *' 01210 print >>n5921reglog,'* 2 N5921_2 15:22:00.00 +05.04.00.00 J2000 *' 01211 print >>n5921reglog,'* Data descriptions: 1 (1 spectral windows and 1 polarization setups) *' 01212 print >>n5921reglog,'* ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs *' 01213 print >>n5921reglog,'* 0 63 LSRK 1412.68608 24.4140625 1550.19688 1413.44902 RR LL *' 01214 print >>n5921reglog,'*********************************' 01215 print "" 01216 print >>n5921reglog,'' 01217 print >>n5921reglog,'********* Export Tests***********' 01218 print >>n5921reglog,'* *' 01219 01220 exportresults = {'UVFITS existence': (uvfitsexists,), 01221 'FITS image existence': (fitsimageexists,)} 01222 if uvfitsexists: 01223 exportresults['UVFITS max'] = (diff_uvfits < 0.05, 01224 '* UVFITS MS max ' + str(fitstest_src), diff_uvfits) 01225 01226 if (fitsimageexists): 01227 exportresults['FITS image max'] = (diff_fitsmax < 0.05, 01228 '* FITS Image max ' + str(fitstest_immax), 01229 diff_fitsmax) 01230 exportresults['FITS image rms'] = (diff_fitsrms < 0.05, 01231 '* FITS Image rms ' + str(fitstest_imrms), 01232 diff_fitsrms) 01233 passedall = reportresults(exportresults) and passedall 01234 01235 print >>n5921reglog,'' 01236 print >>n5921reglog,'********** Regression ***********' 01237 print >>n5921reglog,'* *' 01238 quantresults = {} 01239 quantresults['cal max amplitude'] = (diff_cal < 0.05, 01240 '* Cal max amp ' + str(thistest_cal), diff_cal) 01241 quantresults['src max amplitude'] = (diff_src < 0.05, 01242 '* Src max amp ' + str(thistest_src), diff_src) 01243 quantresults['image max'] = (diff_immax < 0.05, 01244 '* Image max amp ' + str(thistest_immax), diff_immax) 01245 quantresults['image rms'] = (diff_imrms < 0.05, 01246 '* Image rms ' + str(thistest_imrms), diff_imrms) 01247 passedall = reportresults(quantresults) and passedall 01248 01249 if passedall: 01250 regstate=True 01251 print >>n5921reglog,'---' 01252 print >>n5921reglog,'Passed Regression test for NGC5921' 01253 print >>n5921reglog,'---' 01254 print '' 01255 print 'Regression PASSED' 01256 print '' 01257 tstutl.note("Passed Regression test for NGC5921","NORMAL") 01258 else: 01259 regstate=False 01260 print >>n5921reglog,'---' 01261 print >>n5921reglog,'FAILED Regression test for NGC5921' 01262 print >>n5921reglog,'---' 01263 print '' 01264 print 'Regression FAILED' 01265 print '' 01266 tstutl.note("FAILED Regression test for NGC5921","SEVERE") 01267 # Specify what failed here...just saying "FAILED" is aggravating. 01268 for d in (listvisresults, listcalresults, exportresults, quantresults): 01269 listfailures(d) 01270 01271 print >>n5921reglog,'*********************************' 01272 01273 print >>n5921reglog,'' 01274 print >>n5921reglog,'' 01275 print >>n5921reglog,'********* Benchmarking *****************' 01276 print >>n5921reglog,'* *' 01277 print >>n5921reglog,'Total wall clock time was: ', endTime - startTime 01278 print >>n5921reglog,'Total CPU time was: ', endProc - startProc 01279 print >>n5921reglog,'Processing rate MB/s was: ', 35.1/(endTime - startTime) 01280 print >>n5921reglog,'* Breakdown: *' 01281 print >>n5921reglog,'* import time was: '+str(importtime-startTime) 01282 print >>n5921reglog,'* flagautocorr time was: '+str(flagtime-listtime) 01283 print >>n5921reglog,'* setjy time was: '+str(setjytime-flagtime) 01284 print >>n5921reglog,'* bandpass time was: '+str(bptime-setjytime) 01285 print >>n5921reglog,'* gaincal time was: '+str(gaintime-bptime) 01286 print >>n5921reglog,'* listcal time was: '+str(listcaltime-gaintime) 01287 print >>n5921reglog,'* fluxscale time was: '+str(fstime-listcaltime) 01288 print >>n5921reglog,'* applycal time was: '+str(correcttime-fstime) 01289 print >>n5921reglog,'* split-cal time was: '+str(splitcaltime-correcttime) 01290 print >>n5921reglog,'* split-src time was: '+str(splitsrctime-splitcaltime) 01291 print >>n5921reglog,'* contsub time was: '+str(contsubtime-exportuvfitstime) 01292 print >>n5921reglog,'* listvis time was: '+str(listvistime-contsubtime) 01293 print >>n5921reglog,'* clean time was: '+str(cleantime-listvistime) 01294 print >>n5921reglog,'*****************************************' 01295 #print >>n5921reglog,'basho (test cpu) time was: ?? seconds' 01296 01297 print "" 01298 print "Done!" 01299 01300 n5921reglog.close() 01301 01302 #exit() 01303 # 01304 #=====================================================================
1.8.0
