|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 #==================================================================== 00002 # 00003 # Script for simulating and imaging VLA Orion continuum mosaic 00004 # 00005 # Created STM 2008-11-04 from SMA mosaic simulation 00006 # Updated STM 2008-11-06 formatting improvements 00007 # Updated STM 2008-11-19 Patch3 final version 640x640 00008 # Updated STM 2009-06-16 Patch4 800x800 00009 # Updated STM 2009-06-23 Patch4 web download of pickle file 00010 # Updated STM 2009-06-25 minor improvements 00011 # Updated STM 2009-08-24 Release0, new immath syntax, noise 00012 # Updated STM 2009-08-27 pbcor=F on clean 00013 # Updated STM 2009-08-27 noise only option 00014 # Updated STM 2009-10-22 point to Release0 pickle file 00015 # RI 2012-01-15 prefer pickle in DATADIR/regressions/ if exists 00016 # 00017 # Simulates point sources and cleans 00018 # 00019 # Script Notes: 00020 # o Starts with a template MS <templatems> which must be present 00021 # (set <templatems> to point to this) 00022 # o The results are written out as a dictionary in a pickle file 00023 # <scriptprefix>.regression.<datestring>.pickle 00024 # as well as in a text file 00025 # out.<scriptprefix>.<datestring>.log 00026 # (these are not auto-deleted at start of script) 00027 # o If you want to show regression differences from a previous 00028 # run then you can provide a file from a previous run 00029 # <scriptprefix>.pickle 00030 # 00031 #==================================================================== 00032 # IMPORTANT: VERSIONING HERE 00033 #==================================================================== 00034 myscriptvers = '2012-01-15 RI' 00035 00036 import time 00037 import os 00038 import pickle 00039 00040 # needed for existence check 00041 from os import F_OK 00042 00043 scriptprefix='run_orionmos4sim' 00044 prefix='simtest_orionmos4sim' 00045 00046 # Clean up old output files 00047 os.system('rm -rf '+prefix+'*') 00048 00049 #==================================================================== 00050 # IMPORTANT: SET STUFF UP HERE 00051 #==================================================================== 00052 mytmp = '===========================' 00053 mydataset = 'Orion VLA mosaic simulation' 00054 00055 print "========================"+mytmp 00056 print "Simulating and Cleaning "+mydataset 00057 print "========================"+mytmp 00058 print "Script version "+myscriptvers 00059 print "Output will be prefixed with "+prefix 00060 00061 #Start from existing MS 00062 templatems = 'orion_calsplit.ms' 00063 00064 if os.access(templatems,F_OK): 00065 # already in current directory 00066 print " Using "+templatems+" found in current directory" 00067 else: 00068 pathname=os.environ.get('CASAPATH').split()[0] 00069 #datapath=pathname+'/data/regression/ATST3/Orion/' 00070 datapath=pathname+'/data/regression/orionmos4sim/' 00071 msname='orion.ms' 00072 mspath=datapath+msname 00073 # Path to web archive 00074 webms = 'http://casa.nrao.edu/Data/VLA/Orion/'+msname+'.tgz' 00075 # Alternate path at AOC 00076 altms = '/home/ballista/casa/active/data/regression/ATST3/Orion/'+msname 00077 if os.access(msname,F_OK): 00078 # ms in current directory 00079 print " Using "+msname+" found in current directory" 00080 elif os.access(mspath,F_OK): 00081 print '--Copy data to local directory--' 00082 print " Using "+mspath 00083 os.system("cp -r "+mspath+" .") 00084 os.system('chmod -R a+wx '+msname) 00085 elif os.access(altms,F_OK): 00086 print '--Copy data to local directory--' 00087 print " Using "+altms 00088 os.system("cp -r "+altms+" .") 00089 os.system('chmod -R a+wx '+msname) 00090 else: 00091 if os.access(msname+'.tgz',F_OK): 00092 # ms tarball in current directory 00093 print " Using "+msname+".tgz found in current directory" 00094 else: 00095 # try web retrieval 00096 print '--Retrieve data from '+webms 00097 # Use curl (for Mac compatibility) 00098 os.system('curl '+webms+' -f -o '+msname+'.tgz') 00099 # NOTE: could also use wget 00100 #os.system('wget '+webms) 00101 00102 print '--Unpacking tarball ' 00103 os.system('tar xzf '+msname+'.tgz') 00104 if os.access(msname,F_OK): 00105 # should now be in current directory 00106 print " Using "+msname+" found in current directory" 00107 else: 00108 raise IOError," ERROR: "+msname+" not found" 00109 00110 # Starting from orion.ms which has already been calibrated 00111 print '--Split--' 00112 split(vis=msname,outputvis=templatems,datacolumn='corrected') 00113 print " Created "+templatems 00114 00115 #Clean params 00116 myimsize = 800 00117 #myimsize = 640 00118 mycell = "2.0arcsec" 00119 mycen = "J2000 05:35:17.470 -5.23.06.790" 00120 myniter = 400 00121 #mymask = "save.oriontest_mosaic.mask" 00122 mymask = "" 00123 myfield = "2~10" 00124 myspw = "" 00125 myftmachine = "mosaic" 00126 mycyclefactor = 1.5 00127 myminpb = 0.2 00128 00129 #Sim params 00130 myflux = 1.0 00131 mynoise = "" 00132 #mynoise = "0.019Jy" 00133 # EVLA X-band SEFD is 300Jy (1Hz,1sec) so 50MHz,5sec would be 19mJy per vis 00134 00135 print "Will use template MS "+ templatems 00136 print "Will image with imsize %5i and cell %s " % (myimsize,mycell) 00137 print "Will clean "+str(myniter)+" iterations" 00138 if mycen!="": 00139 print "Will use phasecenter "+mycen 00140 if myfield!="": 00141 print "Will image fields "+myfield 00142 if myspw!="": 00143 print "Will use spw "+myspw 00144 if myftmachine!="": 00145 print "Will mosaic using ftmachine "+myftmachine 00146 if mymask!="": 00147 print "Will use clean mask "+mymask 00148 if mynoise!="": 00149 print "Will add noise level of "+mynoise 00150 00151 # Save the parameters used in clean etc. 00152 params = {} 00153 # Version of script 00154 params['version'] = myscriptvers 00155 # User set params 00156 params['user'] = {} 00157 params['user']['templatems'] = templatems 00158 params['user']['myimsize'] = myimsize 00159 params['user']['mycell'] = mycell 00160 params['user']['mycen'] = mycen 00161 params['user']['myniter'] = myniter 00162 params['user']['mymask'] = mymask 00163 params['user']['myfield'] = myfield 00164 params['user']['myspw'] = myspw 00165 params['user']['myftmachine'] = myftmachine 00166 params['user']['mycyclefactor'] = mycyclefactor 00167 params['user']['myminpb'] = myminpb 00168 params['user']['mynoise'] = mynoise 00169 00170 #==================================================================== 00171 00172 # Get this into a direction measure 00173 phc = me.direction(mycen.split()[0],mycen.split()[1],mycen.split()[2]) 00174 # Reference pixel of image 00175 myref = int((myimsize+0.5)/2.0) 00176 dcell = qa.convert(mycell,'arcsec')['value'] 00177 00178 #==================================================================== 00179 # Dictionary with components to simulated 00180 00181 # Comp A is at the peak of the .flux image: 00182 # 05:35:22.024, -05.24.14.789 00183 # Note that only comp A is centered on a pixel by fiat 00184 complist = {} 00185 complist['A'] = {} 00186 complist['A']['rf'] = 'J2000' 00187 complist['A']['RA'] = '05:35:22.024' 00188 complist['A']['Dec'] = '-05.24.14.789' 00189 complist['A']['flux'] = myflux 00190 00191 # Comp B is in center of field 6 (phs center), at 96% of peak .flux 00192 # 05:35:17.47 -05.23.06.79 00193 complist['B'] = {} 00194 complist['B']['rf'] = 'J2000' 00195 complist['B']['RA'] = '05:35:17.47' 00196 complist['B']['Dec'] = '-05.23.06.79' 00197 complist['B']['flux'] = myflux 00198 00199 # Comp C is center of field 10, at 75% of peak .flux 00200 # 05:35:27.52 -05.20.37.52 00201 complist['C'] = {} 00202 complist['C']['rf'] = 'J2000' 00203 complist['C']['RA'] = '05:35:27.52' 00204 complist['C']['Dec'] = '-05.20.37.52' 00205 complist['C']['flux'] = myflux 00206 00207 # Comp D is SSE of field 2, at 33% of peak .flux 00208 # 05:35:08.746 -05.27.19.013 00209 complist['D'] = {} 00210 complist['D']['rf'] = 'J2000' 00211 complist['D']['RA'] = '05:35:08.746' 00212 complist['D']['Dec'] = '-05.27.19.013' 00213 complist['D']['flux'] = myflux 00214 00215 # Which of these to use 00216 compnames = ['A','B','C','D'] 00217 00218 # Box around each component in image 00219 dbox = 10 00220 for comp in compnames: 00221 # get as direction 00222 mydir = me.direction(complist[comp]['rf'],complist[comp]['RA'],complist[comp]['Dec']) 00223 # get offset and angle 00224 off_sep = me.separation(phc,mydir) 00225 off_sec = qa.convert(off_sep,'arcsec')['value'] 00226 off_pa = me.posangle(phc,mydir) 00227 off_par = qa.convert(off_pa,'rad')['value'] 00228 # offsets in pixels 00229 off_x = off_sec*sin(off_par)/dcell 00230 off_y = off_sec*cos(off_par)/dcell 00231 xpix = int( myref + 0.5 - off_x ) 00232 ypix = int( myref + 0.5 + off_y ) 00233 # store these in dictionary 00234 complist[comp]['xpix'] = xpix 00235 complist[comp]['ypix'] = ypix 00236 # now make box for clean and stats 00237 blcx = xpix - dbox 00238 blcy = ypix - dbox 00239 trcx = xpix + dbox 00240 trcy = ypix + dbox 00241 mybox = str(blcx)+","+ str(blcy)+","+ str(trcx)+","+ str(trcy) 00242 complist[comp]['box'] = mybox 00243 myclnbox = [blcx,blcy,trcx,trcy] 00244 complist[comp]['clnbox'] = myclnbox 00245 00246 #==================================================================== 00247 # Start actual processing 00248 #==================================================================== 00249 startTime=time.time() 00250 startProc=time.clock() 00251 00252 #==================================================================== 00253 # Set up component list with 1Jy point source 00254 clfile = prefix + "_sim.cl" 00255 print "Creating componentlist " 00256 00257 for comp in compnames: 00258 mydirection = "J2000 "+complist[comp]['RA']+" "+complist[comp]['Dec'] 00259 compflux = complist[comp]['flux'] 00260 cl.addcomponent(dir=mydirection, flux=compflux, freq='230.0GHz', spectrumtype='constant') 00261 print " Added comp "+comp+" "+str(compflux)+" Jy at "+mydirection+" ("+str(complist[comp]['xpix'])+","+str(complist[comp]['ypix'])+")" 00262 00263 cl.rename(clfile) 00264 mycl = cl.torecord() 00265 cl.close() 00266 00267 print "Created componentlist "+clfile 00268 00269 #Copy a new MS to work from template 00270 msfile = prefix + "_sim.ms" 00271 00272 print "Copying "+templatems+" to "+msfile 00273 00274 os.system('cp -rf '+templatems+' '+msfile) 00275 00276 sm.openfromms(msfile) 00277 sm.setvp(dovp=T, usedefaultvp=T, dosquint=F) 00278 00279 #Add components 00280 sm.predict(complist=clfile) 00281 00282 #Add noise if desired 00283 if mynoise!="": 00284 print 'Adding noise '+mynoise 00285 sm.setvp(dovp=F) 00286 sm.setnoise(mode='simplenoise',simplenoise=mynoise) 00287 sm.corrupt() 00288 00289 sm.close() 00290 00291 print "Completed simulating MS "+msfile 00292 00293 listobs(msfile) 00294 00295 #==================================================================== 00296 #DS9 shape files for viewer 00297 #j2000; text 05:35:07.42 -05:25:36.07 # text={2} color=blue 00298 #j2000; text 05:35:17.42 -05:25:36.79 # text={3} color=blue 00299 #j2000; text 05:35:27.42 -05:25:37.52 # text={4} color=blue 00300 #j2000; text 05:35:27.47 -05:23:07.52 # text={5} color=blue 00301 #j2000; text 05:35:17.47 -05:23:06.79 # text={6} color=blue 00302 #j2000; text 05:35:07.47 -05:23:06.07 # text={7} color=blue 00303 #j2000; text 05:35:07.52 -05:20:36.07 # text={8} color=blue 00304 #j2000; text 05:35:17.52 -05:20:36.80 # text={9} color=blue 00305 #j2000; text 05:35:27.52 -05:20:37.52 # text={10} color=blue 00306 #j2000; text 05:35:32.61 -05:16:07.88 # text={11} color=blue 00307 # 00308 # From listobs: 00309 # Fields: 12 00310 # ID Code Name Right Ascension Declination Epoch 00311 # 0 A 0518+165 05:21:09.89 +16.38.22.04 J2000 00312 # 1 A 0539-057 05:41:38.09 -05.41.49.43 J2000 00313 # 2 ORION1 05:35:07.42 -05.25.36.07 J2000 00314 # 3 ORION2 05:35:17.42 -05.25.36.79 J2000 00315 # 4 ORION3 05:35:27.42 -05.25.37.52 J2000 00316 # 5 ORION4 05:35:27.47 -05.23.07.52 J2000 00317 # 6 ORION5 05:35:17.47 -05.23.06.79 J2000 00318 # 7 ORION6 05:35:07.47 -05.23.06.07 J2000 00319 # 8 ORION7 05:35:07.52 -05.20.36.07 J2000 00320 # 9 ORION8 05:35:17.52 -05.20.36.80 J2000 00321 # 10 ORION9 05:35:27.52 -05.20.37.52 J2000 00322 # 11 ORION10 05:35:32.61 -05.16.07.88 J2000 00323 # Spectral Windows: 00324 # SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 00325 # 0 1 TOPO 8435.1 50000 50000 8435.1 RR LL RL LR 00326 # 1 1 TOPO 8485.1 50000 50000 8485.1 RR LL RL LR 00327 #==================================================================== 00328 00329 sim1time=time.time() 00330 00331 #==================================================================== 00332 # Some details about the mosaic (from old aips++ glish script): 00333 # primary beam at X-band = 5.4', 00334 # mosaic field spacing = 2.5' 00335 # total mosaic size = approx. 9.5' = 570" 00336 # synthesized beam size = 8.4" in D config at 3.6 cm, 8.3 GHz 00337 # cell size = 2" and nx,ny = 300 (600" field size) 00338 # phase center = center field: J2000 05:35:17.470, -005.23.06.790 00339 # NOTE: field 10 is outside of the 9 point primary mosaic (sitting 00340 # on M43 -- but the flux is resolved out so there is no use to 00341 # add it to the mosaic. The script below leaves it out. 00342 #==================================================================== 00343 # Clean 00344 print "Preparing to image data" 00345 00346 # Set up clean boxes 00347 clnbox = [] 00348 for comp in compnames: 00349 clnbox.append(complist[comp]['clnbox']) 00350 00351 print "" 00352 print "Will be cleaning in box(es) "+str(clnbox) 00353 00354 # MFS imaging for continuum 00355 default("clean") 00356 vis = msfile 00357 field = myfield 00358 spw = myspw 00359 selectdata = False 00360 timerange = "" 00361 uvrange = "" 00362 antenna = "" 00363 scan = "" 00364 mode = "mfs" 00365 gain = 0.1 00366 threshold = "0.0mJy" 00367 psfmode = "clark" 00368 imagermode = "mosaic" 00369 ftmachine = myftmachine 00370 mosweight = False 00371 scaletype = "SAULT" 00372 flatnoise = False 00373 multiscale = [] 00374 negcomponent = -1 00375 interactive = F 00376 mask = [] 00377 nchan = 1 00378 start = 0 00379 width = 1 00380 phasecenter = mycen 00381 imsize = myimsize 00382 cell = mycell 00383 restfreq = "" 00384 stokes = "I" 00385 weighting = "natural" 00386 robust = 0.0 00387 uvtaper = False 00388 outertaper = [] 00389 innertaper = [] 00390 modelimage = "" 00391 restoringbeam = [''] 00392 minpb = myminpb 00393 noise = "1.0Jy" 00394 npixels = 0 00395 npercycle = 100 00396 cyclefactor = mycyclefactor 00397 cyclespeedup = -1 00398 00399 # Set phasecenter either on field 1 or using Crystal's 00400 #phasecenter = "1" 00401 print "Will be mosaicing using phasecenter "+phasecenter 00402 00403 print "First make dirty image" 00404 pbcor = T 00405 dirtimag = prefix + ".dirty" 00406 imagename = dirtimag 00407 niter = 0 00408 clean() 00409 print "Output images are called "+imagename 00410 00411 dirty1time=time.time() 00412 00413 if mymask == "": 00414 mask = clnbox 00415 print "Will be cleaning using box(es) "+str(mask) 00416 else: 00417 mask = mymask 00418 print "Will be cleaning using mask "+mask 00419 00420 #Clean image 00421 contimag = prefix + ".clean" 00422 pbcor = F 00423 imagename = contimag 00424 niter = myniter 00425 print "Now actually clean the image using niter = "+str(niter) 00426 clean() 00427 print "Output images are called "+imagename 00428 00429 clean1time=time.time() 00430 00431 # What Crystal's been doing 00432 # 00433 #contimag='g19_d2usb_targets.ms.cont.im' 00434 #clean(vis=vis,imagename=contimag, 00435 # field='',spw='0~12,16~23', 00436 # mode='mfs', 00437 # niter=100,gain=0.1,threshold=0.0, 00438 # psfmode='clark',imagermode='mosaic',scaletype='SAULT', 00439 # mask='g19_d2usb_cont.model.mask', 00440 # interactive=F,imsize=400,cell="0.5arcsec", 00441 # phasecenter="J2000 18h25m56.09 -12d04m28.20", 00442 # pbcor=T,minpb=0.15) 00443 00444 dirimage = dirtimag+'.image' 00445 clnimage = contimag+'.image' 00446 fluximage = contimag+'.flux' 00447 modimage = contimag+'.model' 00448 resimage = contimag+'.residual' 00449 covimage = contimag+'.flux.pbcoverage' 00450 00451 #==================================================================== 00452 # pbcor the model, be careful to mask zeros 00453 00454 pbmodimage = contimag+'.pbcormod' 00455 pbcorimage = contimag+'.pbcorimg' 00456 # Old syntax 00457 #if myflux!=0.0: 00458 # immath(outfile=pbmodimage, 00459 # mode='evalexpr', 00460 # expr="'"+modimage+"'['"+fluximage+"'!=0.0]/'"+fluximage+"'") 00461 #immath(outfile=pbcorimage, 00462 # mode='evalexpr', 00463 # expr="'"+clnimage+"'['"+fluximage+"'!=0.0]/'"+fluximage+"'") 00464 00465 # New syntax (STM 2009-08-24) 00466 if myflux!=0.0: 00467 immath(outfile=pbmodimage, 00468 mode='evalexpr', 00469 imagename=[modimage,fluximage], 00470 expr="IM0[IM1!=0.0]/IM1") 00471 00472 immath(outfile=pbcorimage, 00473 mode='evalexpr', 00474 imagename=[clnimage,fluximage], 00475 expr="IM0[IM1!=0.0]/IM1") 00476 00477 # stats 00478 #==================================================================== 00479 # Some continuum image statistics 00480 00481 # Set up off-src box 00482 # Offset from first component 00483 #offbox = str(myref+45)+','+str(myref+30)+','+str(myref+110)+','+str(myref+100) 00484 doffx = 75 00485 doffy = 65 00486 doffbox = 30 00487 xpix = complist[compnames[0]]['xpix'] 00488 ypix = complist[compnames[0]]['ypix'] 00489 blcx = xpix + doffx - doffbox 00490 blcy = ypix + doffy - doffbox 00491 trcx = xpix + doffx + doffbox 00492 trcy = ypix + doffy + doffbox 00493 offbox = str(blcx)+","+ str(blcy)+","+ str(trcx)+","+ str(trcy) 00494 print "Off-source stats in box "+offbox 00495 00496 # Do dirty image stats 00497 dirstats = imstat(imagename=dirimage,box = '') 00498 dirtyimg_max = dirstats['max'][0] 00499 dirtyimg_rms = dirstats['sigma'][0] 00500 00501 diroffstats = imstat(imagename=dirimage,box=offbox) 00502 dirtyoff_rms = diroffstats['sigma'][0] 00503 00504 # Do clean image 00505 allstats = imstat(imagename=pbcorimage,box = '') 00506 imageall_max = allstats['max'][0] 00507 imageall_rms = allstats['sigma'][0] 00508 00509 offstats = imstat(imagename=pbcorimage,box=offbox) 00510 imageoff_rms = offstats['sigma'][0] 00511 00512 # Do residual image 00513 resstats = imstat(imagename=resimage,box = '') 00514 residall_max = resstats['max'][0] 00515 residall_rms = resstats['sigma'][0] 00516 00517 resoffstats = imstat(imagename=resimage,box=offbox) 00518 residoff_rms = resoffstats['sigma'][0] 00519 00520 if myflux!=0.0: 00521 # Do whole model 00522 modstats = imstat(imagename=modimage,box = '') 00523 modelall_flux = modstats['sum'][0] 00524 00525 pbmodstats = imstat(imagename=pbmodimage,box = '') 00526 pbmodelall_flux = pbmodstats['sum'][0] 00527 00528 # Flux image effective area (STM 2009-06-18) 00529 fluxstats = imstat(imagename=fluximage,box = '') 00530 fluximg_sum = fluxstats['sum'][0] 00531 fluximg_npt = fluxstats['npts'][0] 00532 # Extract pixel area (sq.arcsec.) 00533 pixsec = qa.convert(qa.quantity(mycell),'arcsec') 00534 pixarea = qa.quantity(pixsec)['value']**2 00535 00536 stats1time=time.time() 00537 00538 #==================================================================== 00539 # Some image fitting 00540 # There is a 1Jy source at center of field 1 in image 00541 # Comp Restored image 00542 # Using phasecenter='1' 00543 # A 467,372 at peak between 2,3,4 18:25:55.116 -12.04.32.500 00544 # B 400,400 at 75% center field 1 18:25:57.400 -12.04.18.50 00545 # C 360,462 at 40% N of field 0 18:25:58.764 -12.03.47.504 00546 # D 511,420 at 20% NW of field 3 18:25:53.607 -12.04.08.582 00547 00548 # Using phasecenter="J2000 18h25m56.09 -12d04m28.20" 00549 # Note that only comp A is centered on a pixel by fiat 00550 # A 429,391 at peak between 2,3,4 18:25:55.101 -12.04.32.700 00551 # B 362,420 at 75% center field 1 18:25:57.400 -12.04.18.50 00552 # C 321,481 at 40% N of field 0 18:25:58.764 -12.03.47.504 00553 # D 473,440 at 20% NW of field 3 18:25:53.607 -12.04.08.582 00554 00555 # Use 21x21 boxes centered on above coordinates 00556 # Calculate: 00557 # (a) Peak in restored image (Jy/bm) 00558 # (b) Flux in imfit of restored image (Jy) 00559 # (c) Flux in box in sim model (Jy) 00560 # 00561 00562 totalflux = 0.0 00563 # Box around each component in image 00564 for comp in compnames: 00565 compflux = complist[comp]['flux'] 00566 totalflux = totalflux + compflux 00567 mybox = complist[comp]['box'] 00568 print "Component "+comp+" stats in box "+mybox 00569 # Dirty image 00570 xstat_dirimg = imstat(imagename=dirimage,box=mybox) 00571 # Clean image 00572 xstat_img = imstat(imagename=pbcorimage,box=mybox) 00573 if myflux!=0.0: 00574 # Clean model 00575 xstat_mod = imstat(imagename=modimage,box=mybox) 00576 xstat_pbmod = imstat(imagename=pbmodimage,box=mybox) 00577 # Residual image 00578 xstat_res = imstat(imagename=resimage,box=mybox) 00579 # Save in complist 00580 complist[comp]['dirtystat'] = xstat_dirimg 00581 complist[comp]['imstat'] = xstat_img 00582 if myflux!=0.0: 00583 complist[comp]['modstat'] = xstat_mod 00584 complist[comp]['pbmodstat'] = xstat_pbmod 00585 complist[comp]['resstat'] = xstat_res 00586 00587 #==================================================================== 00588 # Done 00589 #==================================================================== 00590 endProc=time.clock() 00591 endTime=time.time() 00592 00593 import datetime 00594 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00595 00596 myvers = casalog.version() 00597 myuser = os.getenv('USER') 00598 myhost = str( os.getenv('HOST') ) 00599 mycwd = os.getcwd() 00600 myos = os.uname() 00601 mypath = os.environ.get('CASAPATH') 00602 00603 walltime = (endTime - startTime) 00604 cputime = (endProc - startProc) 00605 00606 # Print to terminal, and also save most things to a logfile 00607 outfile = 'out.'+prefix+'.'+datestring+'.log' 00608 logfile = open(outfile,'w') 00609 00610 def lprint(msg, lfile): 00611 """ 00612 Prints msg to both stdout and lfile. 00613 """ 00614 print msg 00615 print >> logfile, msg 00616 00617 lprint('Running '+myvers+' on host '+myhost, logfile) 00618 lprint(' at '+datestring, logfile) 00619 lprint(' using '+mypath, logfile) 00620 00621 # 00622 ########################################################################## 00623 # Write some more stuff to logfile 00624 lprint('', logfile) 00625 lprint('---', logfile) 00626 lprint('Script version: '+params['version'], logfile) 00627 lprint('User set parameters used in execution:', logfile) 00628 lprint('---', logfile) 00629 00630 for keys in params['user'].keys(): 00631 lprint(' %s = %s ' % ( keys, params['user'][keys] ), logfile) 00632 00633 lprint('', logfile) 00634 00635 # 00636 ########################################################################## 00637 # 00638 # Save info in regression dictionary 00639 new_regression = {} 00640 00641 new_regression['date'] = datestring 00642 new_regression['version'] = myvers 00643 new_regression['user'] = myuser 00644 new_regression['host'] = myhost 00645 new_regression['cwd'] = mycwd 00646 new_regression['os'] = myos 00647 new_regression['aipspath'] = mypath 00648 00649 new_regression['dataset'] = mydataset 00650 00651 # Store the parameters used in clean etc. 00652 new_regression['params'] = params 00653 00654 # Store the results of this regression 00655 results = {} 00656 00657 # First the overall results 00658 op = 'divf' 00659 tol = 0.01 00660 results['clean_dirtyimg_max'] = {} 00661 results['clean_dirtyimg_max']['name'] = 'Dirty image (all) max' 00662 results['clean_dirtyimg_max']['value'] = dirtyimg_max 00663 results['clean_dirtyimg_max']['op'] = op 00664 results['clean_dirtyimg_max']['tol'] = tol 00665 00666 results['clean_dirtyimg_rms'] = {} 00667 results['clean_dirtyimg_rms']['name'] = 'Dirty image (all) rms' 00668 results['clean_dirtyimg_rms']['value'] = dirtyimg_rms 00669 results['clean_dirtyimg_rms']['op'] = op 00670 results['clean_dirtyimg_rms']['tol'] = tol 00671 00672 results['clean_dirtyoff_rms'] = {} 00673 results['clean_dirtyoff_rms']['name'] = 'Dirty image (off-src) rms' 00674 results['clean_dirtyoff_rms']['value'] = dirtyoff_rms 00675 results['clean_dirtyoff_rms']['op'] = op 00676 results['clean_dirtyoff_rms']['tol'] = tol 00677 00678 results['clean_imageall_max'] = {} 00679 results['clean_imageall_max']['name'] = 'Clean image (all) max' 00680 results['clean_imageall_max']['value'] = imageall_max 00681 results['clean_imageall_max']['op'] = op 00682 results['clean_imageall_max']['tol'] = tol 00683 00684 results['clean_imageall_rms'] = {} 00685 results['clean_imageall_rms']['name'] = 'Clean image (all) rms' 00686 results['clean_imageall_rms']['value'] = imageall_rms 00687 results['clean_imageall_rms']['op'] = op 00688 results['clean_imageall_rms']['tol'] = tol 00689 00690 results['clean_imageoff_rms'] = {} 00691 results['clean_imageoff_rms']['name'] = 'Clean image (off-src) rms' 00692 results['clean_imageoff_rms']['value'] = imageoff_rms 00693 results['clean_imageoff_rms']['op'] = op 00694 results['clean_imageoff_rms']['tol'] = tol 00695 00696 results['clean_residall_max'] = {} 00697 results['clean_residall_max']['name'] = 'Residual image (all) max' 00698 results['clean_residall_max']['value'] = residall_max 00699 results['clean_residall_max']['op'] = op 00700 results['clean_residall_max']['tol'] = tol 00701 00702 results['clean_residall_rms'] = {} 00703 results['clean_residall_rms']['name'] = 'Residual image (all) rms' 00704 results['clean_residall_rms']['value'] = residall_rms 00705 results['clean_residall_rms']['op'] = op 00706 results['clean_residall_rms']['tol'] = tol 00707 00708 results['clean_residoff_rms'] = {} 00709 results['clean_residoff_rms']['name'] = 'Residual image (off-src) rms' 00710 results['clean_residoff_rms']['value'] = residoff_rms 00711 results['clean_residoff_rms']['op'] = op 00712 results['clean_residoff_rms']['tol'] = tol 00713 00714 if myflux!=0.0: 00715 results['clean_modelall_flux'] = {} 00716 results['clean_modelall_flux']['name'] = 'Model image (all) flux' 00717 results['clean_modelall_flux']['value'] = modelall_flux 00718 results['clean_modelall_flux']['op'] = op 00719 results['clean_modelall_flux']['tol'] = tol 00720 00721 results['clean_pbmodelall_flux'] = {} 00722 results['clean_pbmodelall_flux']['name'] = 'PBCOR Model image (all) flux' 00723 results['clean_pbmodelall_flux']['value'] = pbmodelall_flux 00724 results['clean_pbmodelall_flux']['op'] = op 00725 results['clean_pbmodelall_flux']['tol'] = tol 00726 00727 results['sim_comp_flux'] = {} 00728 results['sim_comp_flux']['name'] = 'Sim Model total flux' 00729 results['sim_comp_flux']['value'] = totalflux 00730 results['sim_comp_flux']['op'] = op 00731 results['sim_comp_flux']['tol'] = tol 00732 00733 results['flux_area'] = {} 00734 results['flux_area']['name'] = 'Eff. mosaic area (sq.arcmin.)' 00735 results['flux_area']['value'] = fluximg_sum*pixarea/3600.0 00736 results['flux_area']['op'] = op 00737 results['flux_area']['tol'] = tol 00738 00739 # Done filling results 00740 new_regression['results'] = results 00741 00742 # Now the component results 00743 new_regression['compnames'] = compnames 00744 #Store raw statistics 00745 new_regression['compstats'] = complist 00746 00747 compresults = {} 00748 for comp in compnames: 00749 compresults[comp] = {} 00750 00751 xstat_dirimg = complist[comp]['dirtystat'] 00752 xstat_img = complist[comp]['imstat'] 00753 if myflux!=0.0: 00754 xstat_mod = complist[comp]['modstat'] 00755 xstat_pbmod = complist[comp]['pbmodstat'] 00756 xstat_res = complist[comp]['resstat'] 00757 00758 compresults[comp]['box'] = complist[comp]['box'] 00759 00760 op = 'divf' 00761 tol = 0.08 00762 compresults[comp]['dirty_image_max'] = {} 00763 compresults[comp]['dirty_image_max']['name'] = 'Dirty image max' 00764 compresults[comp]['dirty_image_max']['value'] = xstat_dirimg['max'][0] 00765 compresults[comp]['dirty_image_max']['op'] = op 00766 compresults[comp]['dirty_image_max']['tol'] = tol 00767 00768 op = 'diff' 00769 tol = 0.5 00770 compresults[comp]['dirty_image_maxposx'] = {} 00771 compresults[comp]['dirty_image_maxposx']['name'] = 'Dirty image max pos x' 00772 compresults[comp]['dirty_image_maxposx']['value'] = xstat_dirimg['maxpos'][0] 00773 compresults[comp]['dirty_image_maxposx']['op'] = op 00774 compresults[comp]['dirty_image_maxposx']['tol'] = tol 00775 00776 compresults[comp]['dirty_image_maxposy'] = {} 00777 compresults[comp]['dirty_image_maxposy']['name'] = 'Dirty image max pos y' 00778 compresults[comp]['dirty_image_maxposy']['value'] = xstat_dirimg['maxpos'][1] 00779 compresults[comp]['dirty_image_maxposy']['op'] = op 00780 compresults[comp]['dirty_image_maxposy']['tol'] = tol 00781 00782 op = 'divf' 00783 tol = 0.08 00784 compresults[comp]['clean_image_max'] = {} 00785 compresults[comp]['clean_image_max']['name'] = 'Restored image max' 00786 compresults[comp]['clean_image_max']['value'] = xstat_img['max'][0] 00787 compresults[comp]['clean_image_max']['op'] = op 00788 compresults[comp]['clean_image_max']['tol'] = tol 00789 00790 op = 'diff' 00791 tol = 0.5 00792 compresults[comp]['clean_image_maxposx'] = {} 00793 compresults[comp]['clean_image_maxposx']['name'] = 'Restored image max pos x' 00794 compresults[comp]['clean_image_maxposx']['value'] = xstat_img['maxpos'][0] 00795 compresults[comp]['clean_image_maxposx']['op'] = op 00796 compresults[comp]['clean_image_maxposx']['tol'] = tol 00797 00798 compresults[comp]['clean_image_maxposy'] = {} 00799 compresults[comp]['clean_image_maxposy']['name'] = 'Restored image max pos y' 00800 compresults[comp]['clean_image_maxposy']['value'] = xstat_img['maxpos'][1] 00801 compresults[comp]['clean_image_maxposy']['op'] = op 00802 compresults[comp]['clean_image_maxposy']['tol'] = tol 00803 00804 op = 'divf' 00805 tol = 0.08 00806 compresults[comp]['residual_image_max'] = {} 00807 compresults[comp]['residual_image_max']['name'] = 'Residual image max' 00808 compresults[comp]['residual_image_max']['value'] = xstat_res['max'][0] 00809 compresults[comp]['residual_image_max']['op'] = op 00810 compresults[comp]['residual_image_max']['tol'] = tol 00811 00812 compresults[comp]['residual_image_rms'] = {} 00813 compresults[comp]['residual_image_rms']['name'] = 'Residual image rms' 00814 compresults[comp]['residual_image_rms']['value'] = xstat_res['sigma'][0] 00815 compresults[comp]['residual_image_rms']['op'] = op 00816 compresults[comp]['residual_image_rms']['tol'] = tol 00817 00818 if myflux!=0.0: 00819 compresults[comp]['model_image_sum'] = {} 00820 compresults[comp]['model_image_sum']['name'] = 'Model image flux' 00821 compresults[comp]['model_image_sum']['value'] = xstat_mod['sum'][0] 00822 compresults[comp]['model_image_sum']['op'] = op 00823 compresults[comp]['model_image_sum']['tol'] = tol 00824 00825 compresults[comp]['pbmodel_image_sum'] = {} 00826 compresults[comp]['pbmodel_image_sum']['name'] = 'PBCOR Model image flux' 00827 compresults[comp]['pbmodel_image_sum']['value'] = xstat_pbmod['sum'][0] 00828 compresults[comp]['pbmodel_image_sum']['op'] = op 00829 compresults[comp]['pbmodel_image_sum']['tol'] = tol 00830 00831 compresults[comp]['sim_comp_flux'] = {} 00832 compresults[comp]['sim_comp_flux']['name'] = 'Simulated comp flux' 00833 compresults[comp]['sim_comp_flux']['value'] = complist[comp]['flux'] 00834 compresults[comp]['sim_comp_flux']['op'] = op 00835 compresults[comp]['sim_comp_flux']['tol'] = tol 00836 00837 op = 'diff' 00838 tol = 0.5 00839 compresults[comp]['sim_comp_posx'] = {} 00840 compresults[comp]['sim_comp_posx']['name'] = 'Simulated comp pos x' 00841 compresults[comp]['sim_comp_posx']['value'] = complist[comp]['xpix'] 00842 compresults[comp]['sim_comp_posx']['op'] = op 00843 compresults[comp]['sim_comp_posx']['tol'] = tol 00844 00845 compresults[comp]['sim_comp_posy'] = {} 00846 compresults[comp]['sim_comp_posy']['name'] = 'Simulated comp pos y' 00847 compresults[comp]['sim_comp_posy']['value'] = complist[comp]['ypix'] 00848 compresults[comp]['sim_comp_posy']['op'] = op 00849 compresults[comp]['sim_comp_posy']['tol'] = tol 00850 00851 new_regression['compresults'] = compresults 00852 00853 # Dataset size info 00854 datasize_raw = 110.0 # MB 00855 datasize_ms = 110.0 # MB 00856 new_regression['datasize'] = {} 00857 new_regression['datasize']['raw'] = datasize_raw 00858 new_regression['datasize']['ms'] = datasize_ms 00859 00860 # Save timing to regression dictionary 00861 new_regression['timing'] = {} 00862 00863 total = {} 00864 total['wall'] = (endTime - startTime) 00865 total['cpu'] = (endProc - startProc) 00866 total['rate_raw'] = (datasize_raw/(endTime - startTime)) 00867 total['rate_ms'] = (datasize_ms/(endTime - startTime)) 00868 00869 new_regression['timing']['total'] = total 00870 00871 nstages = 4 00872 new_regression['timing']['nstages'] = nstages 00873 00874 stages = {} 00875 stages[0] = ['simulate',(sim1time-startTime)] 00876 stages[1] = ['invert',(dirty1time-sim1time)] 00877 stages[2] = ['clean',(clean1time-dirty1time)] 00878 stages[3] = ['stats',(endTime-clean1time)] 00879 00880 new_regression['timing']['stages'] = stages 00881 00882 # 00883 ########################################################################## 00884 # See if a previous pickfile can be found 00885 # 00886 # Try and load previous results from regression file 00887 # 00888 regression = {} 00889 regressfile = scriptprefix + '.pickle' 00890 prev_results = {} 00891 00892 repodir=os.getenv("CASAPATH") 00893 repodir=repodir.split()[0] 00894 regressdirfile= repodir+ "/data/regression/orionmos4sim/"+scriptprefix + '.pickle' 00895 00896 # Path to web archive (latest version) 00897 #webfile = 'http://casa.nrao.edu/Doc/Scripts/'+scriptprefix+'.pickle' 00898 #webfile = 'http://casa.nrao.edu/Patch4/Doc/Scripts/'+scriptprefix+'.pickle' 00899 webfile = 'http://casa.nrao.edu/Release0/Doc/Scripts/'+scriptprefix+'.pickle' 00900 00901 if os.access(regressfile,F_OK): 00902 # pickle file in current directory 00903 print " Found "+regressfile+" in current directory" 00904 else: 00905 if os.access(regressdirfile,F_OK): 00906 # pickle file in local copy of regressions dir 00907 print " Found "+regressdirfile+" in your local data repo" 00908 regressfile=regressdirfile 00909 else: 00910 print "Trying web download of "+webfile 00911 os.system('curl '+webfile+' -f -o '+regressfile) 00912 # NOTE: could also use wget 00913 # os.system('wget '+webfile) 00914 00915 # Now try to access this file 00916 try: 00917 fr = open(regressfile,'r') 00918 except: 00919 print "No previous regression results file "+regressfile 00920 regression['exist'] = False 00921 else: 00922 u = pickle.Unpickler(fr) 00923 regression = u.load() 00924 fr.close() 00925 print "Regression results filled from "+regressfile 00926 print "Regression from version "+regression['version']+" on "+regression['date'] 00927 regression['exist'] = True 00928 00929 prev_results = regression['results'] 00930 00931 prev_compresults = regression['compresults'] 00932 00933 # 00934 ########################################################################## 00935 # Now print regression results if a previous pickfile can be found 00936 if myflux!=0.0: 00937 resultlist = ['clean_dirtyimg_max','clean_dirtyimg_rms', 00938 'clean_imageall_max','clean_imageall_rms','clean_imageoff_rms', 00939 'clean_residall_max','clean_residall_rms','clean_residoff_rms', 00940 'clean_modelall_flux','clean_pbmodelall_flux', 00941 'sim_comp_flux','flux_area'] 00942 00943 compreslist = ['dirty_image_max','dirty_image_maxposx','dirty_image_maxposy', 00944 'clean_image_max','clean_image_maxposx','clean_image_maxposy', 00945 'residual_image_max','residual_image_rms', 00946 'model_image_sum','pbmodel_image_sum', 00947 'sim_comp_flux','sim_comp_posx','sim_comp_posy'] 00948 else: 00949 resultlist = ['clean_dirtyimg_max','clean_dirtyimg_rms', 00950 'clean_imageall_max','clean_imageall_rms','clean_imageoff_rms', 00951 'clean_residall_max','clean_residall_rms','clean_residoff_rms', 00952 'flux_area'] 00953 00954 compreslist = ['dirty_image_max','dirty_image_maxposx','dirty_image_maxposy', 00955 'clean_image_max','clean_image_maxposx','clean_image_maxposy', 00956 'residual_image_max','residual_image_rms'] 00957 00958 testresults = {} 00959 testcompresults = {} 00960 00961 if regression['exist']: 00962 lprint('---', logfile) 00963 lprint('Regression versus previous values:', logfile) 00964 lprint('---', logfile) 00965 lprint(" Regression results filled from "+regressfile, logfile) 00966 lprint(" Regression from version " + regression['version'] + " on " + 00967 regression['date'], logfile) 00968 lprint(" Regression platform " + regression['host'] + " using " + 00969 regression['aipspath'], logfile) 00970 final_status = 'Passed' 00971 else: 00972 # Just print results 00973 lprint('---', logfile) 00974 lprint('Results of imaging:', logfile) 00975 lprint('---', logfile) 00976 lprint(" No previous regression file", logfile) 00977 00978 # Do the component results 00979 for comp in compnames: 00980 lprint("", logfile) 00981 lprint(" Component " + comp + " :", logfile) 00982 00983 if regression['exist']: 00984 # Set up storage for regression results 00985 testcompresults[comp] = {} 00986 00987 for keys in compreslist: 00988 res = compresults[comp][keys] 00989 new_val = res['value'] 00990 00991 if regression['exist']: 00992 testcompresults[comp][keys] = {} 00993 testres = {} 00994 testres['value'] = new_val 00995 testres['op'] = res['op'] 00996 testres['tol'] = res['tol'] 00997 00998 if prev_compresults.has_key(comp): 00999 # This is a known comp 01000 if prev_compresults[comp].has_key(keys): 01001 # This is a known regression key 01002 prev = prev_compresults[comp][keys] 01003 prev_val = prev['value'] 01004 if res['op'] == 'divf': 01005 new_diff = (new_val - prev_val)/prev_val 01006 else: 01007 new_diff = new_val - prev_val 01008 01009 if abs(new_diff) > res['tol']: 01010 new_status = 'Failed' 01011 else: 01012 new_status = 'Passed' 01013 01014 testres['prev'] = prev_val 01015 testres['diff'] = new_diff 01016 testres['status'] = new_status 01017 testres['test'] = 'Last' 01018 01019 else: 01020 # Unknown regression key 01021 testres['prev'] = 0.0 01022 testres['diff'] = 1.0 01023 testres['status'] = 'Missed' 01024 testres['test'] = 'none' 01025 01026 else: 01027 # Unknown regression comp 01028 testres['prev'] = 0.0 01029 testres['diff'] = 1.0 01030 testres['status'] = 'Unknown' 01031 testres['test'] = 'none' 01032 01033 # Save results 01034 testcompresults[comp][keys] = testres 01035 01036 if testres['status']=='Failed': 01037 final_status = 'Failed' 01038 01039 # Print results 01040 lprint('--%30s : %12.6f was %12.6f %4s %12.6f (%6s) %s ' % (res['name'], 01041 testres['value'], 01042 testres['prev'], 01043 testres['op'], 01044 testres['diff'], 01045 testres['status'], 01046 testres['test']), 01047 logfile) 01048 else: 01049 # Just print results 01050 lprint('--%30s : %12.6f ' % ( res['name'], res['value'] ), logfile) 01051 01052 # Do the overall results 01053 lprint("", logfile) 01054 lprint(" Whole image :", logfile) 01055 for keys in resultlist: 01056 res = results[keys] 01057 new_val = res['value'] 01058 if regression['exist']: 01059 # Set up storage for regression results 01060 testresults[keys] = {} 01061 testres = {} 01062 testres['value'] = new_val 01063 testres['op'] = res['op'] 01064 testres['tol'] = res['tol'] 01065 01066 # Compute regression results 01067 if prev_results.has_key(keys): 01068 # This is a known regression key 01069 prev = prev_results[keys] 01070 prev_val = prev['value'] 01071 if res['op'] == 'divf': 01072 new_diff = (new_val - prev_val)/prev_val 01073 else: 01074 new_diff = new_val - prev_val 01075 01076 if abs(new_diff)>res['tol']: 01077 new_status = 'Failed' 01078 else: 01079 new_status = 'Passed' 01080 01081 testres['prev'] = prev_val 01082 testres['diff'] = new_diff 01083 testres['status'] = new_status 01084 testres['test'] = 'Last' 01085 else: 01086 # Unknown regression key 01087 testres['prev'] = 0.0 01088 testres['diff'] = 1.0 01089 testres['status'] = 'Missed' 01090 testres['test'] = 'none' 01091 01092 # Save results 01093 testresults[keys] = testres 01094 01095 # Print results 01096 lprint('--%30s : %12.6f was %12.6f %4s %12.6f (%6s) %s ' % (res['name'], 01097 testres['value'], 01098 testres['prev'], 01099 testres['op'], 01100 testres['diff'], 01101 testres['status'], 01102 testres['test'] ), 01103 logfile) 01104 if testres['status']=='Failed': 01105 final_status = 'Failed' 01106 01107 else: 01108 # Just print the current results 01109 lprint('--%30s : %12.6f ' % ( res['name'], res['value'] ), logfile) 01110 01111 # Done 01112 01113 if regression['exist']: 01114 # Final regression status 01115 if (final_status == 'Passed'): 01116 regstate=True 01117 print >>logfile,'---' 01118 lprint('Passed Regression test for ' + mydataset, logfile) 01119 print >>logfile,'---' 01120 print '' 01121 print 'Regression PASSED' 01122 print '' 01123 else: 01124 regstate=False 01125 lprint('----FAILED Regression test for ' + mydataset, logfile) 01126 01127 # Write these regression test results to the dictionary 01128 new_regression['testresults'] = testresults 01129 new_regression['testcompresults'] = testcompresults 01130 else: 01131 regstate=False 01132 lprint('----Unable to complete regression test for '+ mydataset, logfile) 01133 01134 # 01135 ########################################################################## 01136 # Final stats and timing 01137 01138 lprint('', logfile) 01139 lprint('********* Benchmarking *************************', logfile) 01140 lprint('* *', logfile) 01141 lprint('Total wall clock time was: %10.3f ' % total['wall'], logfile) 01142 lprint('Total CPU time was: %10.3f ' % total['cpu'], logfile) 01143 lprint('Raw processing rate MB/s was: %8.1f ' % total['rate_raw'], logfile) 01144 lprint('MS processing rate MB/s was: %8.1f ' % total['rate_ms'], logfile) 01145 lprint('', logfile) 01146 lprint('* Breakdown: *', logfile) 01147 01148 for i in range(nstages): 01149 lprint('* %16s * time was: %10.3f ' % tuple(stages[i]), logfile) 01150 01151 lprint('************************************************', logfile) 01152 if regression['exist']: 01153 lprint('Previous wall time was %10.3f sec on %s ' % 01154 (regression['timing']['total']['wall'], regression['host'] ), logfile) 01155 lprint('Previous CPU time was %10.3f sec on %s ' % 01156 (regression['timing']['total']['cpu'], regression['host'] ), logfile) 01157 01158 # 01159 ########################################################################## 01160 # 01161 # Save regression results as dictionary using Pickle 01162 # 01163 pickfile = 'out.'+prefix + '.regression.'+datestring+'.pickle' 01164 f = open(pickfile,'w') 01165 p = pickle.Pickler(f) 01166 p.dump(new_regression) 01167 f.close() 01168 01169 print "" 01170 print "Regression result dictionary saved in "+pickfile 01171 print "" 01172 print "Use Pickle to retrieve these" 01173 print "" 01174 01175 # e.g. 01176 # f = open(pickfile) 01177 # u = pickle.Unpickler(f) 01178 # clnmodel = u.load() 01179 # polmodel = u.load() 01180 # f.close() 01181 01182 lprint("", logfile) 01183 lprint("Done with " + mydataset, logfile) 01184 01185 logfile.close() 01186 print "Results are in "+outfile 01187 01188 # 01189 ##########################################################################
1.8.0
