|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################### 00002 # # 00003 # Regression/Benchmarking Script for Orion # 00004 # Single Dish + Synthesis # 00005 ############################################### 00006 00007 import time 00008 import os 00009 import shutil 00010 00011 os.system('rm -rf orion_t* gbt_gau.im') 00012 00013 pathname=os.environ.get('CASAPATH') 00014 datapath=os.environ.get('CASAPATH').split()[0]+'/data/regression/ATST3/Orion/' 00015 00016 print '--Copy data to local directory--' 00017 mspath='cp -r '+datapath+'orion.ms .' 00018 os.system(mspath) 00019 os.system('chmod -R a+wx orion.ms') 00020 00021 startTime = time.time() 00022 startProc = time.clock() 00023 00024 00025 print '--Feather--' 00026 # Starting from: 00027 # VLA Orion mosaic image : orion_vlamem.im 00028 # GBT OTF image : orion.gbt.im 00029 feather('orion_tfeather.im',datapath+'orion_vlamem.im',datapath+'orion.gbt.im') 00030 feathertime = time.time() 00031 #GBT: Max:5.129806e+01 Flux:2.425065e+02 Jy rms:1.277546e+01 00032 #VLA: Max:8.340111e-01 Flux:1.891523e+02 Jy rms:1.099514e-01 00033 00034 print '--Feather - create synth image--' 00035 # Starting from: 00036 # VLA Orion mosaic image : orion.ms 00037 # GBT OTF image : orion.gbt.im 00038 # Some details about the mosaic: 00039 # primary beam at X-band = 5.4', 00040 # mosaic field spacing = 2.5' 00041 # total mosaic size = approx. 9.5' = 570" 00042 # synthesized beam size = 8.4" in D config at 3.6 cm, 8.3 GHz 00043 # cell size = 2" and nx,ny = 300 (600" field size) 00044 # phase center = center field: J2000 05:35:17.470, -005.23.06.790 00045 # NOTE: field 10 is outside of the 9 point primary mosaic (sitting 00046 # on M43 -- but the flux is resolved out so there is no use to 00047 # add it to the mosaic. The script below leaves it out. 00048 shutil.rmtree('orion.task.model', True) 00049 default('mosaic') 00050 mosaic('orion.ms', 00051 'orion.task', 00052 'mfs', 00053 'entropy', 00054 niter=35, 00055 sigma='4mJy', 00056 targetflux='240Jy', 00057 mask=datapath+'orion.mask6', 00058 field=[2, 3, 4, 5, 6, 7, 8, 9,10], 00059 spw=[0, 1], 00060 stokes='I', 00061 cell=[2, 2], 00062 imsize=[300, 300], 00063 weighting='briggs', 00064 mosweight=True, 00065 rmode='norm', 00066 robust=-1, 00067 cyclefactor=4, 00068 cyclespeedup=500, 00069 prior='', 00070 phasecenter=6, 00071 ftmachine='ft', 00072 minpb=0.1, 00073 scaletype='PBCOR') 00074 00075 feather('orion_tfeather2.im','orion.task.image',datapath+'orion.gbt.im') 00076 #GBT: Max:5.129806e+01 Flux:2.425065e+02 Jy rms:1.277546e+01 00077 #VLA: Max:8.340111e-01 Flux:1.891523e+02 Jy rms:1.099514e-01 00078 feathersynthtime = time.time() 00079 00080 print '--Single Dish as Model (multi-scale)--' 00081 ## Starting from: 00082 ## VLA calibrated visibilities: orion.ms 00083 ## GBT OTF cube: orion.gbt.im 00084 shutil.rmtree('orion_tsdms.model', True) 00085 shutil.rmtree('orion_tsdms.mask', True) 00086 default('clean') 00087 clean(vis='orion.ms', 00088 imagename='orion_tsdms', 00089 field='2~10', 00090 spw='0,1', 00091 mode='mfs', 00092 niter=10000, 00093 gain=0.2, 00094 threshold='10.0mJy', 00095 imagermode='mosaic', 00096 mosweight=True, 00097 ftmachine='ft', 00098 cyclefactor=4, 00099 cyclespeedup=-1, 00100 multiscale=[0,2,8,16,32,64], 00101 negcomponent=-1, 00102 mask=datapath+'orion.mask6', 00103 modelimage=datapath+'orion.gbt.im', 00104 imsize=[300,300], 00105 cell=['2.0arcsec','2.0arcsec'], 00106 phasecenter=6, 00107 stokes='I', 00108 weighting='briggs', 00109 robust=-1.0, 00110 pbcor=True, 00111 minpb=0.1) 00112 sdmodelmstime = time.time() 00113 ###combo: Max:1.195286e+00 Flux:2.873779e+02 Jy rms:9.069330e-02 00114 ###GBT: Max:5.129806e+01 Flux:2.425065e+02 Jy rms:1.277546e+01 00115 ###VLA: Max:8.340111e-01 Flux:1.891523e+02 Jy rms:1.099514e-01 00116 ## 00117 print '--Single Dish as Model (MEM)--' 00118 ### Starting from: 00119 ### VLA calibrated visibilities: orion.ms 00120 ### GBT OTF cube: orion.gbt.im 00121 shutil.rmtree('orion_tsdmem.model', True) 00122 default('mosaic') 00123 mosaic('orion.ms', 'orion_tsdmem', 'mfs', 'entropy', niter=3, sigma='4mJy', 00124 targetflux='240Jy', mask=datapath+'orion.mask6', 00125 field=[2, 3, 4, 5, 6, 7, 8, 9, 10], spw=[0, 1], stokes='I', 00126 cell=[2, 2], imsize=[300, 300], 00127 weighting='briggs', mosweight=True, rmode='norm', robust=-1, 00128 cyclefactor=4, cyclespeedup=500, 00129 phasecenter=6, 00130 modelimage='orion_tsdmem', sdimage=datapath+'orion.gbt.im', 00131 ftmachine='ft', 00132 prior='orion_tsdmem', 00133 minpb=0.1, scaletype='PBCOR') 00134 sdmodelmemtime=time.time() 00135 ###GBT: Max:5.129806e+01 Flux:2.425065e+02 Jy rms:1.277546e+01 00136 ###VLA: Max:8.340111e-01 Flux:1.891523e+02 Jy rms:1.099514e-01 00137 ## 00138 ###### 00139 #### Create synthesis (BIMA) data cube, deconvolve single dish cube 00140 #### DO joint deconvolution 00141 ###### 00142 def joint_deconvolve(datapath): 00143 print '--Joint deconvolution --' 00144 00145 #Regrid GBT image onto synth imaging coordinates 00146 ia.open('orion_tsdms.image') 00147 csys = ia.coordsys() 00148 ia.close() 00149 ia.open(datapath+'orion.gbt.im') 00150 ib=ia.regrid(outfile='orion_tgbt_regrid.im',shape=[300,300,1,1], axes=[0,1], 00151 csys=csys.torecord(),overwrite=True) 00152 ib.setcoordsys(csys.torecord()) 00153 ia.close() 00154 ib.done() 00155 00156 #Deconvolve GBT image 00157 # Sigh. dc.open will warn about the lack of a PSF, but I can't seem to 00158 # define a PSF before calling dc.open. 00159 dc.open('orion_tgbt_regrid.im', psf='') 00160 #make gaussian for PSF (best guess for GBT beam based on beamsize 00161 #report in GBT image) 00162 dc.makegaussian('gbt_gau.im',bmaj='55arcsec',bmin='55arcsec', 00163 bpa='0deg',normalize=False) 00164 dc.close() 00165 os.system("rm -rf orion_tjoint3") 00166 dc.open('orion_tgbt_regrid.im',psf='gbt_gau.im') 00167 dc.setscales(scalemethod='uservector',uservector=[30.,100.,200.]) 00168 dc.clean(algorithm='fullmsclean',model='orion_tjoint3',niter=50, 00169 gain=0.3,mask='',threshold='0.5Jy') 00170 dc.close() 00171 #default('clean') 00172 ia.open('orion_tjoint3') 00173 ia.calc(pixels='orion_tjoint3*"'+datapath+'orion.mask6"') 00174 ia.close() 00175 00176 im.open('orion.ms') 00177 im.selectvis(field=[2,3,4,5,6,7,8,9,10],spw=[0,1]) 00178 im.defineimage(nx=300,cellx='2arcsec',phasecenter=6,spw=[0,1]) 00179 im.setvp(dovp=True) 00180 #im.setscales(scalemethod='uservector',uservector=[0,3,10,30,100]) 00181 im.setscales(scalemethod='uservector',uservector=[0,3,10,30]) 00182 ###if clean component for large scale goes negative continue to use 00183 ##that scale 00184 im.setoptions(ftmachine='ft') 00185 im.setmfcontrol(stoplargenegatives=-1, cyclefactor=2.0, cyclespeedup=-1) 00186 im.weight(type='briggs',rmode='norm',robust=-1,mosaic=True) 00187 im.clean(algorithm='mfmultiscale', model='orion_tjoint3', 00188 image='orion_tjoint3.image', gain=0.3, niter=10000, 00189 mask=datapath+'orion.mask6', threshold='4mJy') 00190 im.close() 00191 return time.time() 00192 00193 jointmemtime = joint_deconvolve(datapath) 00194 00195 endProc = time.clock() 00196 endTime = time.time() 00197 00198 # Regression 00199 import datetime 00200 datestring = datetime.datetime.isoformat(datetime.datetime.today()) 00201 outfile = 'orion.' + datestring + '.log' 00202 logfile = open(outfile, 'w') 00203 00204 test_results = {} 00205 for k, fn in (('Feather 1', 'orion_tfeather.im'), 00206 ('Feather 2', 'orion_tfeather2.im'), 00207 ('SD Model (MS)', 'orion_tsdms.image'), 00208 ('SD Model (MEM)', 'orion_tsdmem.image'), 00209 ('Joint Deconvolution', 'orion_tjoint3.image')): 00210 if ia.open(fn): 00211 test_results[k] = ia.statistics() 00212 ia.close() 00213 else: 00214 print >>logfile, "Could not open", fn, "for reading!" 00215 00216 print >>logfile,'' 00217 print >>logfile,'********** Data Summary *********' 00218 print >>logfile,' GBT image ' 00219 print >>logfile,'Image name : orion.gbt.im' 00220 print >>logfile,'Object name :' 00221 print >>logfile,'Image type : PagedImage' 00222 print >>logfile,'Image quantity : Intensity' 00223 print >>logfile,'Pixel mask(s) : mask0' 00224 print >>logfile,'Region(s) : None' 00225 print >>logfile,'Image units : Jy/beam' 00226 print >>logfile,'Restoring Beam : 98.0547 arcsec, 98.0547 arcsec, 79.25 deg' 00227 print >>logfile,'Direction reference : J2000' 00228 print >>logfile,'Spectral reference : LSRK' 00229 print >>logfile,'Velocity type : RADIO' 00230 print >>logfile,'Rest frequency : 1 Hz' 00231 print >>logfile,'Telescope : GBT' 00232 print >>logfile,'Observer : UNKNOWN' 00233 print >>logfile,'Date observation : UNKNOWN' 00234 print >>logfile,'Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units' 00235 print >>logfile,'------------------------------------------------------------------------------------------------' 00236 print >>logfile,'1 1 Direction Right Ascension SIN 300 300 05:35:17.470 151.00 -2.000000e+00 arcsec' 00237 print >>logfile,'2 1 Direction Declination SIN 300 300 -05.23.06.790 151.00 2.000000e+00 arcsec' 00238 print >>logfile,'3 2 Stokes Stokes 1 1 I' 00239 print >>logfile,'4 3 Spectral Frequency 1 1 8.4351e+09 1.00 6.050000e+07 Hz' 00240 print >>logfile,' Velocity -inf 1.00 nan km/s' 00241 print >>logfile,' Observer: unavailable Project: DSTST' 00242 print >>logfile,'Observation: VLA(26 antennas)' 00243 print >>logfile,'Data records: 1093716 Total integration time = 8545 seconds' 00244 print >>logfile,' Observed from 11:15:48 to 13:38:13' 00245 print >>logfile,'Fields: 12' 00246 print >>logfile,' ID Name Right Ascension Declination Epoch' 00247 print >>logfile,' 0 0518+165 05:21:09.89 +16.38.22.04 J2000' 00248 print >>logfile,' 1 0539-057 05:41:38.09 -05.41.49.43 J2000' 00249 print >>logfile,' 2 ORION1 05:35:07.42 -05.25.36.07 J2000' 00250 print >>logfile,' 3 ORION2 05:35:17.42 -05.25.36.79 J2000' 00251 print >>logfile,' 4 ORION3 05:35:27.42 -05.25.37.52 J2000' 00252 print >>logfile,' 5 ORION4 05:35:27.47 -05.23.07.52 J2000' 00253 print >>logfile,' 6 ORION5 05:35:17.47 -05.23.06.79 J2000' 00254 print >>logfile,' 7 ORION6 05:35:07.47 -05.23.06.07 J2000' 00255 print >>logfile,' 8 ORION7 05:35:07.52 -05.20.36.07 J2000' 00256 print >>logfile,' 9 ORION8 05:35:17.52 -05.20.36.80 J2000' 00257 print >>logfile,' 10 ORION9 05:35:27.52 -05.20.37.52 J2000' 00258 print >>logfile,' 11 ORION10 05:35:32.61 -05.16.07.88 J2000' 00259 print >>logfile,'Spectral Windows: (2 unique spectral windows and 1 unique polarization setups)' 00260 print >>logfile,' SpwID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs' 00261 print >>logfile,' 0 1 TOPO 8435.1 50000 50000 8435.1 RR RL LR LL' 00262 print >>logfile,' 1 1 TOPO 8485.1 50000 50000 8485.1 RR RL LR LL' 00263 print >>logfile,'*********************************' 00264 print >>logfile,'' 00265 print >>logfile,'********** Regression ***********' 00266 print >>logfile,'* *' 00267 00268 # Test name Stat type Expected Label irregularities 00269 test_descs = (('Feather 1', 'max', 0.780, ' '), 00270 ('Feather 2', 'max', 0.978, ' '), 00271 ('SD Model (MS)', 'max', 1.14), 00272 ('SD Model (MEM)', 'max', 0.906), 00273 ('Joint Deconvolution', 'max', 1.10, '', 'Joint Decon1'), # 1.014 00274 ('Feather 1', 'flux', 242.506, ' '), 00275 ('Feather 2', 'flux', 242.506, ' '), 00276 ('SD Model (MS)', 'flux', 347, ' ', 'SD Model (MS)', 'SD Model MS'), 00277 ('SD Model (MEM)', 'flux', 286, '', 'SD Model (MEM)', 'Joint Deconvolution'), 00278 ('Joint Deconvolution', 'flux', 226, '', 'Joint Decon2')) # 360.468 00279 00280 def log_test_result(test_results, testdesc, logfile): 00281 """Append testdesc to logfile and return whether or not the test was 00282 successful.""" 00283 result = test_results[testdesc[0]][testdesc[1]][0] 00284 reldiff = abs(1.0 - result / testdesc[2]) 00285 if reldiff < 0.05: 00286 print >>logfile, '* Passed', 00287 print '* Alright', 00288 retval = True 00289 else: 00290 print >>logfile, '! FAILED', 00291 print '#####FAILED' 00292 retval = False 00293 00294 # RR 4/18/2009: I think this complication might stem from a bug in the 00295 # original version of the script, but since it is a regression script I 00296 # am hesitant to change the output. 00297 if len(testdesc) > 5: 00298 title1 = testdesc[5] 00299 else: 00300 title1 = testdesc[0] 00301 if len(testdesc) > 4: 00302 title2 = testdesc[4] 00303 else: 00304 title2 = testdesc[0] 00305 title2 += ':' 00306 if testdesc[1] == 'max': 00307 title1 += ' image max' 00308 title2 += ' Image max' 00309 else: 00310 title1 += ' ' + testdesc[1] 00311 title2 += ' ' + testdesc[1].title() 00312 if len(testdesc) > 3: 00313 title2 += testdesc[3] 00314 00315 print >>logfile, title1, 'test' 00316 print title1, 'test' 00317 print >>logfile, '*-- ' + title2 + str(result) + ',' + str(testdesc[2]) 00318 print '*-- ' + title2 + str(result) + ',' + str(testdesc[2]) 00319 return retval 00320 00321 regstate = True 00322 for td in test_descs: 00323 regstate &= log_test_result(test_results, td, logfile) 00324 00325 if regstate: 00326 print >>logfile, '---' 00327 print >>logfile, 'Passed Regression test for Orion' 00328 print >>logfile, '---' 00329 print '' 00330 print 'Regression PASSED' 00331 print '' 00332 else: 00333 print >>logfile, '----FAILED Regression test for Orion' 00334 print '' 00335 print 'Regression FAILED' 00336 print '' 00337 00338 print >>logfile, '*********************************' 00339 00340 print >>logfile,'' 00341 print >>logfile,'' 00342 print >>logfile,'********* Benchmarking *****************' 00343 print >>logfile,'* *' 00344 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00345 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00346 print >>logfile,'Processing rate MB/s was: '+str(760./(endTime - startTime)) 00347 print >>logfile,'* Breakdown: *' 00348 print >>logfile,'* feather time was: '+str(feathertime-startTime) 00349 print >>logfile,'* feathersynth time was: '+str(feathersynthtime-feathertime) 00350 print >>logfile,'* sdmodel (MS) time was: '+str(sdmodelmstime-feathersynthtime) 00351 print >>logfile,'* sdmodel(MEM) time was: '+str(sdmodelmemtime-sdmodelmstime) 00352 print >>logfile,'* joint decon time was: '+str(jointmemtime-sdmodelmemtime) 00353 print >>logfile,'*****************************************' 00354 # 00355 logfile.close()
1.8.0
