|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################ 00002 # test clean spw and channelization in various ways 00003 # This version will test only the following tests defined in tests[{}] 00004 # avel,nvel,chan,freq,frm,fr.5. 00005 # In order to run everything, uncomment the other tests in tests[{}] (line 147), 00006 # but it will take many hours. 00007 00008 00009 import os, time 00010 00011 # options you can set in local variables: 00012 # pyonly: just run through setChannelization, not full clean 00013 # analonly: print stats about previously run data and make plot. 00014 00015 # e.g. 00016 # CASA> pyonly=True 00017 # CASA> execfile("cleanchan_regression.py") 00018 00019 00020 # root name of files produced: 00021 rt="cln_w3oh" 00022 pyonly=False 00023 00024 l=locals() 00025 if not l.has_key("analonly"): 00026 analonly=False 00027 if not analonly: 00028 os.system("rm -rf "+rt+"*") 00029 00030 startTime = time.time() 00031 startProc = time.clock() 00032 00033 print '--Running clean chan/spw test--' 00034 00035 import datetime 00036 datestring = datetime.datetime.isoformat(datetime.datetime.today()) 00037 outfile = rt + datestring + '.log' 00038 logfile = open(outfile, 'w') 00039 print 'Writing output to ' + outfile + "\n" 00040 00041 l=locals() 00042 if not l.has_key("repodir"): 00043 repodir=os.getenv("CASAPATH").split(' ')[0] 00044 00045 print 'I think the data repository is at '+repodir 00046 datadir=repodir+"/data/regression/cvel/input/" 00047 00048 # just do python part, not full on clean 00049 if not l.has_key("pyonly"): 00050 pyonly=False 00051 00052 # get ms 00053 if not os.path.exists(rt+".ms"): 00054 importuvfits(fitsfile=datadir+'W3OH_MC.UVFITS', vis=rt+'.ms') 00055 00056 00057 # get stats of ms 00058 tb.open(rt+'.ms/SPECTRAL_WINDOW') 00059 cf = tb.getvarcol('CHAN_FREQ') 00060 ch0_ms = cf['r1'][0][0] / 1e3 00061 chn_ms = cf['r1'][510][0] / 1e3 00062 del_ms = (cf['r1'][1][0] - 1e3*ch0_ms) /1e3 00063 del2_ms = (1e3*chn_ms - cf['r1'][509][0]) /1e3 00064 cw = tb.getvarcol('CHAN_WIDTH') 00065 wid_ms = cw['r1'][0][0] /1e3 00066 nch_ms = tb.getcol("NUM_CHAN")[0] 00067 reffreq_ms = tb.getcol("REF_FREQUENCY")[0] / 1e3 00068 tb.close() 00069 00070 # get line peak in ms: 00071 tb.open(rt+".ms") 00072 dat=tb.getcol("DATA") 00073 tb.done() 00074 spec=pl.zeros(nch_ms) 00075 for i in range(nch_ms): 00076 spec[i]=pl.mean(dat[0,i,]) 00077 00078 # the line is around ch260 00079 x=pl.array(range(21))+250 00080 y=spec[x] 00081 fitchan_ms = pl.mean((x+1)*y)/pl.mean(y) 00082 x0=floor(fitchan_ms) 00083 xf=fitchan_ms-x0 00084 fit_ms = cf['r1'][x0][0] + xf * (cf['r1'][x0+1][0]-cf['r1'][x0][0]) 00085 fit_ms=fit_ms / 1e3 00086 00087 del dat 00088 # manually set the line center - the above doesn't work well because of the blended lines 00089 fit_ms = 1665658.566 00090 00091 # PLOT 00092 #f=cf['r1'] 00093 #f=f.flatten() 00094 # 00095 #pl.ion() 00096 #pl.clf() 00097 #pl.plot(f,spec/4.6,label='ms') 00098 #pl.plot(f,spec/4.6,'mo',label='ms') 00099 #pl.xlim([1665640000,1665675000]) 00100 00101 00102 00103 00104 # start accumulating big array of stats 00105 if pyonly: 00106 stats=[(nch_ms,"%fkHz" % ch0_ms,"%fkHz" % del_ms)] 00107 else: 00108 stats=[(ch0_ms,del_ms,wid_ms,del2_ms,nch_ms,chn_ms,fit_ms)] 00109 00110 # and the name of the test associated with each row 00111 lname=['ms'] 00112 sname=['ms'] 00113 00114 00115 00116 00117 # function to analyze images 00118 def imstats(image): 00119 ia.open(image) 00120 ch0 = (ia.toworld([n/2,n/2,0,0],'n')['numeric'][3])/ 1e3 00121 nch = ia.shape()[3] 00122 chn = (ia.toworld([n/2,n/2,0,nch-1],'n')['numeric'][3])/ 1e3 00123 del1 = (ia.toworld([n/2,n/2,0,1],'n')['numeric'][3] - 1e3*ch0) /1e3 00124 wid = (ia.summary()['incr'][3]) /1e3 00125 del2 = (1e3*chn - ia.toworld([n/2,n/2,0,nch-2],'n')['numeric'][3]) /1e3 00126 try: 00127 mylist = [] 00128 for i in range(len(ia.shape())): 00129 mylist.append(0) 00130 mylist.append(0) 00131 t = tuple(mylist) 00132 center = ia.fitprofile(ngauss=1,poly=1)['gs']['center'].item(t) 00133 csys = ia.coordsys() 00134 fit = csys.velocitytofrequency(center)/1e3 00135 except: 00136 fit=0. 00137 return (ch0,del1,wid,del2,nch,chn,fit) 00138 00139 # function to swap first/last for vel results, for new convention of vel means 00140 # _increasing_ vel by default 00141 def swapvel(stats): 00142 return (stats[5], -stats[1], -stats[2], -stats[3], stats[4], stats[0], stats[6]) 00143 00144 00145 00146 00147 # start/end velocities in ms: 00148 f0=reffreq_ms 00149 v1_ms=(f0-ch0_ms)/f0*299792.458 00150 v0_ms=(f0-chn_ms)/f0*299792.458 00151 incr = -del_ms/f0*299792.458 00152 w4="%fkm/s" % incr 00153 00154 # tests of various modes: 00155 tests=[ 00156 {'name':'avel', 'mode':'velocity', 'desc':'velo default'}, 00157 # {'name':'vel', 'mode':'velocity', 'spw':'0:200~299', 'desc':'velo spw=0:200~299'}, 00158 # {'name':'vel3', 'mode':'velocity', 'spw':'0:200~299','sta':50, 'desc':'velo spw=0:200~299 sta=50'}, 00159 # {'name':'vel4', 'mode':'velocity', 'spw':'0:200~299','sta':250, 'desc':'velo spw=0:200~299 sta=250'}, 00160 # {'name':'vel2', 'mode':'velocity', 'spw':'0:210~250,0:280~290', 'desc':'velo spw=0:210~250,0:280~290'}, 00161 # {'name':'cha3', 'mode':'channel', 'spw':'0:210~250,0:280~290', 'desc':'chan spw=0:210~250,0:280~290'}, 00162 # {'name':'sve1', 'mode':'velocity', 'sta':"%fkm/s" % v0_ms, 'desc':'velo start=%5.2fkm/s' % v0_ms}, 00163 # {'name':'sve2', 'mode':'velocity', 'sta':"%fkm/s" % (v0_ms-0.01), 'desc':'velo start=%5.2fkm/s' % (v0_ms-0.01)}, 00164 # {'name':'sve3', 'mode':'velocity', 'sta':"%fkm/s" % (v0_ms-0.01), 'desc':'velo start=%5.2fkm/s linear' % (v0_ms-0.01), 'int':'linear'}, 00165 # {'name':'sve4', 'mode':'velocity', 'sta':"%fkm/s" % (v0_ms-0.1), 'desc':'velo start=%5.2fkm/s' % (v0_ms-0.1)}, 00166 # {'name':'sve5', 'mode':'velocity', 'sta':"%fkm/s" % -15.105, 'desc':'velo start=%5.2fkm/s' % -15.105}, 00167 {'name':'nvel', 'mode':'velocity', 'sta':"%fkm/s" % v1_ms,'wid':w4,'desc':'velo start=%5.2fkm/s wid=%5.2fkm/s' % (v1_ms,incr)}, 00168 # {'name': 'fel', 'mode':'velocity', 'veltype':'optical', 'desc':'felo default'}, 00169 {'name':'chan', 'mode':'channel', 'desc':'chan default'}, 00170 # {'name':'cha1', 'mode':'channel', 'sta':200,'wid':10, 'desc':'chan 200w10'}, 00171 # {'name':'cha2', 'mode':'channel', 'sta':200,'nchan':99, 'desc':'chan 200-299'}, 00172 # {'name':'cha4', 'mode':'channel', 'shift':wid_ms/3, 'desc':'chan start=%5.3fkHz' % (wid_ms/3)}, 00173 # {'name':'chfr', 'mode':'channel', 'wid':2.5, 'desc':'chan wid=2.5'}, 00174 # {'name':'chfm', 'mode':'channel', 'wid':-1, 'desc':'chan wid=-1'}, # segfaults 00175 # {'name':'frch', 'mode':'frequency','wid':2.5, 'desc':'freq wid=2.5'}, 00176 {'name':'frm', 'mode':'frequency','wid':"-%fkHz" % wid_ms,'sta':'%fkHz' % chn_ms, 'desc':'freq wid=-%6.1fKHz' % wid_ms}, 00177 {'name':'freq', 'mode':'frequency', 'desc':'freq default'}, 00178 # {'name':'fr.3', 'mode':'frequency','shift':wid_ms/3, 'desc':'freq start=%5.3fkHz' % (wid_ms/3)}, 00179 {'name':'fr.5', 'mode':'frequency','shift':wid_ms/2, 'desc':'freq start=%5.3fkHz' % (wid_ms/2)}, 00180 # {'name':'fr.7', 'mode':'frequency','shift':wid_ms/3*2, 'desc':'freq start=%5.3fkHz' % (wid_ms*2/3)}, 00181 # {'name':'fr.3l','mode':'frequency','shift':wid_ms/3,'int':'linear','desc':'freq start=%5.3fkHz lin' % (wid_ms/3)} 00182 ] 00183 00184 00185 # a subset for the plot 00186 #toplot=['avel','fel','chan','freq','fr.5'] 00187 toplot=[] 00188 00189 if pyonly: 00190 from cleanhelper import * 00191 imCln=imtool() 00192 vis=rt+'.ms' 00193 imset=cleanhelper(imCln, vis, False, casalog) 00194 outframe='' 00195 else: 00196 # size of image to make 00197 n=64 00198 00199 00200 #debug: 00201 #tests=[tests[18]] 00202 00203 00204 j=0 00205 for te in tests: 00206 if te.has_key('spw'): 00207 spw=te['spw'] 00208 else: 00209 spw='' 00210 if te.has_key('nchan'): 00211 nchan=te['nchan'] 00212 else: 00213 nchan=-1 00214 if te.has_key('int'): 00215 interp=te['int'] 00216 else: 00217 interp='nearest' 00218 if te.has_key('wid'): 00219 wid=te['wid'] 00220 else: 00221 wid='' 00222 # should not be required but is in old versions of cleanhelper: 00223 # if te['mode']=='channel': 00224 # wid=1 00225 00226 if te.has_key('veltype'): 00227 vtype=te['veltype'] 00228 else: 00229 vtype='radio' 00230 00231 if te.has_key('sta'): 00232 sta=te['sta'] 00233 elif te.has_key('shift'): 00234 start=ch0_ms + te['shift'] 00235 sta="%f kHz" % start 00236 else: 00237 sta='' 00238 # should not be required but is in old versions of cleanhelper: 00239 # if te['mode']=='channel': 00240 # sta=0 00241 00242 00243 if pyonly: 00244 try: 00245 print te['mode'],spw,"",nchan,sta,wid,outframe,vtype,str(reffreq_ms)+'kHz' 00246 # test the new version 00247 st = imset.setChannelizeDefault(te['mode'],spw,"",nchan,sta,wid,outframe,vtype,"", str(reffreq_ms)+'kHz') 00248 #st = imset.setChannelization(te['mode'],spw,"",nchan,sta,wid,outframe,vtype,str(reffreq_ms)+'kHz') 00249 # (localnchan, localstart, localwidth)=imset.setChannelization(mode,spw,field,nchan,start,width,outframe,veltype,restfreq) 00250 except: 00251 st = (-1,"","") 00252 else: 00253 try: 00254 if not analonly: 00255 print 'Running clean for name=%s, mode=%s'%(te['name'],te['mode']) 00256 clean(vis=rt+'.ms', 00257 imagename=rt+'_'+te['name'], 00258 cell="6arcsec",imsize=[n,n], 00259 spw=spw, 00260 imagermode='', 00261 mode=te['mode'],veltype=vtype, 00262 start=sta,width=wid,nchan=nchan, 00263 interpolation=interp, 00264 niter=100,threshold='1mJy', 00265 restfreq=str(reffreq_ms)+'kHz') 00266 if te['mode']=='velocity': 00267 st=swapvel(imstats(rt+'_'+te['name']+'.image')) 00268 else: 00269 st=imstats(rt+'_'+te['name']+'.image') 00270 if te.has_key('shift'): 00271 # shift back: 00272 shift=te['shift'] 00273 st=(st[0]-shift, st[1],st[2],st[3],st[4], st[5]-shift, st[6]-shift) 00274 00275 if toplot.count(te['name'])>0: 00276 ia.open(rt+'_'+te['name']+'.image') 00277 foo=ia.getchunk(blc=[25,25,0,0],trc=[40,40,0,510],dropdeg=True,axes=[0,1]) 00278 c=ia.coordsys() 00279 n=len(foo) 00280 f=pl.array(range(n)) 00281 for i in range(n): 00282 f[i]=c.toworld([32,32,0,i])['numeric'][3] 00283 ia.done() 00284 #f=f+100*j 00285 pl.plot(f,foo,label=te['desc']) 00286 j=j+1 00287 except: 00288 st=(reffreq_ms,0,0,0,-1,reffreq_ms,reffreq_ms) 00289 stats.append(st) 00290 lname.append(te['desc']) 00291 sname.append(te['name']) 00292 00293 00294 00295 00296 # max name len 00297 lth=0 00298 for i in range(len(lname)): 00299 lth0=len(lname[i]) 00300 if lth0>lth: 00301 lth=lth0 00302 fmt="%-"+str(lth)+"s" 00303 00304 # header line 00305 if pyonly: 00306 yo = fmt+" nchan chan0 width " 00307 else: 00308 yo = fmt+" chan0 (kHz) ch1-ch0 width chn-(n-1) n lastchan fit peak" 00309 00310 00311 # print "regular" stats: 00312 # print yo % " " 00313 print >> logfile, yo % " " 00314 for i in range(len(lname)): 00315 if pyonly: 00316 #print fmt % sname[i], "%4i %20s %20s" % stats[i] 00317 print >> logfile, fmt % lname[i], "%4i %20s %20s" % stats[i] 00318 else: 00319 print >> logfile, fmt % lname[i], "%15.7f %11.7f %10.7f %10.7f %4i %16.7f %16.7f" % stats[i] 00320 00321 00322 #from matplotlib.font_manager import fontManager, FontProperties 00323 #font= FontProperties(size='x-small') 00324 #pl.legend(prop=font,loc=2) 00325 #pl.xlim([1665632000, 1665675000]) 00326 #pl.ylim([-5,100]) 00327 #pl.savefig( rt + datestring + ".png") 00328 00329 00330 # print stats w/chans as pixel fractions relative to MS: 00331 00332 print yo % " " 00333 if pyonly: 00334 del imCln 00335 regstate=True 00336 # print ms values: 00337 print fmt % lname[0], "%4i %20s %20s" % stats[0] 00338 # ms goes from ch0_ms to chn_ms in freq. 00339 for i in range(1,len(lname)): 00340 if stats[i][0]>0: 00341 if tests[i-1]['mode']=='velocity': 00342 if tests[i-1].has_key('veltype'): 00343 vtype=tests[i-1]['veltype'] 00344 else: 00345 vtype='radio' 00346 v0=qa.convert(stats[i][1],'km/s')['value'] 00347 v1=v0+qa.convert(stats[i][2],'km/s')['value'] 00348 if vtype=='radio': 00349 f0 = reffreq_ms * ( 1 - v0/299792.458) 00350 f1 = reffreq_ms * ( 1 - v1/299792.458) 00351 else: 00352 f0 = reffreq_ms / ( 1 + v0/299792.458) 00353 f1 = reffreq_ms / ( 1 + v1/299792.458) 00354 w=f1-f0 00355 elif tests[i-1]['mode']=='channel': 00356 if type(stats[i][1])==str: 00357 f0=qa.convert(stats[i][1],'kHz')['value'] 00358 w=qa.convert(stats[i][2],'kHz')['value'] 00359 else: 00360 f0=ch0_ms + wid_ms*stats[i][1] 00361 w=wid_ms*stats[i][2] 00362 else: 00363 if type(stats[i][1])==str: 00364 f0 = qa.convert(stats[i][1],'kHz')['value'] 00365 w = qa.convert(stats[i][2],'kHz')['value'] 00366 else: 00367 f0=ch0_ms + wid_ms*stats[i][1] 00368 w=wid_ms*stats[i][2] 00369 print fmt % lname[i], "%4i ch %15f %16fkHz" % (stats[i][0], (f0-ch0_ms)/wid_ms, w) 00370 print >> logfile, fmt % lname[i], "%4i ch %15f %16fkHz" % (stats[i][0], (f0-ch0_ms)/wid_ms, w) 00371 else: 00372 print fmt % lname[i], " FAIL" 00373 print >> logfile, fmt % lname[i], " FAIL" 00374 else: 00375 print fmt % lname[0], "%15.7f %10.7f %10.7f %10.7f %4i %16.7f %16.7f" % stats[0] 00376 print >> logfile, fmt % lname[0], "%15.7f %10.7f %10.7f %10.7f %4i %16.7f %16.7f" % stats[0] 00377 for i in range(1,len(lname)): 00378 foo = pl.array(stats[i]) 00379 if stats[i][4]>0: 00380 foo[0] = (foo[0]-ch0_ms)/wid_ms 00381 foo[5] = (foo[5]-ch0_ms)/wid_ms 00382 foo[6] = (foo[6]-fit_ms)/wid_ms 00383 print fmt % lname[i], "%15.7f %10.7f %10.7f %10.7f %4i %16.7f %16.7f" % (foo[0],foo[1],foo[2],foo[3],foo[4],foo[5],foo[6]) 00384 print >> logfile, fmt % lname[i], "%15.7f %10.7f %10.7f %10.7f %4i %16.7f %16.7f" % (foo[0],foo[1],foo[2],foo[3],foo[4],foo[5],foo[6]) 00385 else: 00386 print fmt % lname[i], " FAIL" 00387 print >> logfile, fmt % lname[i], " FAIL" 00388 00389 00390 # regress 00391 regstate=True 00392 00393 # test everything against ms assumed to be first stat 00394 # tolerance 00395 tol=0.005 00396 00397 tests=['ch0','ch0-ch1','width','chn-ch(n-1)','nchan','chn','fit'] 00398 for run in range(1,len(lname)-1): 00399 for te in range(len(tests)): 00400 adiff=abs(stats[run][te]-stats[0][te])/stats[0][te] 00401 if tests[te]=='width' or tests[te]=='ch0-ch1' or tests[te]=='chn-ch(n-1)': 00402 if sname[run]=='nvel' or sname[run]=='frm': 00403 adiff=abs(-stats[run][te]-stats[0][te])/stats[0][te] 00404 if adiff < tol: 00405 print >> logfile, "* Passed %10s %-10s test, got %-11.5g expected %-11.5g" % (sname[run],tests[te],stats[run][te],stats[0][te]) 00406 else: 00407 print >> logfile, "* FAILED %10s %-10s test, got %-11.5g expected %-11.5g " % (sname[run],tests[te],stats[run][te],stats[0][te]) 00408 regstate = False 00409 00410 00411 00412 00413 print >> logfile,'---' 00414 if regstate: 00415 print >> logfile, 'Passed', 00416 print '' 00417 print 'Regression PASSED' 00418 print '' 00419 else: 00420 print >> logfile, 'FAILED', 00421 print '' 00422 print 'Regression FAILED' 00423 print '' 00424 print >> logfile, 'regression test for clean chan/spw' 00425 print >>logfile,'---' 00426 print >>logfile,'*********************************' 00427 00428 endTime = time.time() 00429 endProc = time.clock() 00430 00431 print >>logfile,'' 00432 print >>logfile,'********** Benchmarking **************' 00433 print >>logfile,'' 00434 print >>logfile,'Total wall clock time was: %8.3f s.' % (endTime - startTime) 00435 print >>logfile,'Total CPU time was: %8.3f s.' % (endProc - startProc) 00436 print >>logfile,'Wall processing rate was: %8.3f MB/s.' % (17896.0 / 00437 (endTime - startTime)) 00438 00439 00440 logfile.close() 00441 00442 00443 print '--Finished clean chan/spw, written to '+outfile+' --'
1.8.0
