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