|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################################################# 00002 # $Id:$ 00003 # Test Name: # 00004 # Regression Test Script for the cvel task # 00005 # # 00006 # Rationale for Inclusion: # 00007 # The task cvel needs to be exercised and compared to clean # 00008 # # 00009 # Features tested: # 00010 # 1) does cvel run without raising exceptions for input frame TOPO # 00011 # and all possible output frames? # 00012 # 2) can clean process the cvel output? # 00013 # 3) does cvel+clean produce compatible results to clean-only? # 00014 # (channel flux values, channel world coordinates) # 00015 # # 00016 # Input data: # 00017 # one dataset, one scan of a VLA observation provided by Crystal Brogan # 00018 # # 00019 ############################################################################# 00020 00021 # copy module needed to create dictionaries to store the results 00022 import copy 00023 00024 myname = 'cvel_regression' 00025 00026 # default dataset name 00027 dataset_name_orig = 'ANTEN_sort_hann_for_cvel_reg-thinned.ms' 00028 00029 # get the dataset name from the wrapper if possible 00030 mydict = locals() 00031 if mydict.has_key("dataset_name"): 00032 dataset_name_orig = mydict["dataset_name"] 00033 00034 def isnear(a,b,p): 00035 print " ", a, b 00036 if(a==b): 00037 print " exactly equal" 00038 return True 00039 dev = abs(a-b) 00040 print " deviation = ", dev 00041 if(dev<=p): 00042 return True 00043 return False 00044 00045 def isrnear(a,b,p): 00046 print " ", a, b 00047 if(a==b): 00048 print " exactly equal" 00049 return True 00050 rdev = abs(a-b)/abs(max(a,b)) 00051 print " relative deviation = ", rdev 00052 if(rdev<=p): 00053 return True 00054 return False 00055 00056 00057 # start frequencies of the grids for the different output reference frames 00058 00059 freqmodestart = { 'TOPO': '2.21078e+10Hz', 00060 'LSRK': '2.21067e+10Hz', 00061 'LSRD': '2.21066e+10Hz', 00062 'BARY': '2.21065e+10Hz', 00063 'GALACTO': '2.21181e+10Hz', 00064 'LGROUP': '2.2125e+10Hz', 00065 'CMB': '2.20804e+10Hz' 00066 } 00067 00068 # hanning smooth switches for the different output frames 00069 00070 dohanning = { 'TOPO': False, 00071 'LSRK': True, 00072 'LSRD': False, 00073 'BARY': False, 00074 'GALACTO': False, 00075 'LGROUP': False, 00076 'CMB': False 00077 } 00078 00079 #channel widths 00080 00081 freqmodewidth = { 'TOPO': '1.953125e+05Hz', 00082 'LSRK': '1.953027e+05Hz', 00083 'LSRD': '1.953017e+05Hz', 00084 'BARY': '1.953005e+05Hz', 00085 'GALACTO': '1.954028e+05Hz', 00086 'LGROUP': '1.954639e+05Hz', 00087 'CMB': '1.950702e+05Hz' 00088 } 00089 00090 # channel to test with imstats 00091 00092 peakchan = '28' 00093 otherchan1 = '5' 00094 otherchan2 = '25' 00095 otherchan3 = '30' 00096 otherchan4 = '45' 00097 00098 testregion = '135,127,135,127' 00099 00100 # storage for results 00101 imstats = { 'TOPO': 0,'LSRK': 0, 'LSRD': 0, 'BARY': 0, 'GALACTO': 0, 'LGROUP': 0, 'CMB': 0 } 00102 mode_imstats = { peakchan: copy.deepcopy(imstats), 00103 otherchan1: copy.deepcopy(imstats), 00104 otherchan2: copy.deepcopy(imstats), 00105 otherchan3: copy.deepcopy(imstats), 00106 otherchan4: copy.deepcopy(imstats) } 00107 cvel_imstats = { 'frequency': copy.deepcopy(mode_imstats), 00108 'radio velocity': copy.deepcopy(mode_imstats), 00109 'optical velocity': copy.deepcopy(mode_imstats), 00110 'channel': copy.deepcopy(mode_imstats) } 00111 cleanonly_imstats = copy.deepcopy(cvel_imstats) 00112 00113 00114 # First we want to test the regridding and spectral frame transformation only. 00115 # So we combine the spws in the input dataset beforehand. 00116 # Also: clean needs write access to the input MS, so we need a local copy anyway. 00117 00118 dataset_name = dataset_name_orig 00119 os.system('rm -rf input.ms input2.ms') 00120 os.system('cp -RL '+dataset_name_orig+' input.ms') 00121 os.system('chmod -R u+w input.ms') 00122 os.system('cp -RL '+dataset_name_orig+' input2.ms') 00123 os.system('chmod -R u+w input2.ms') 00124 clean_inputvis_local_copy = 'input.ms' 00125 00126 clean_inputvis_local_copy2 = 'input2.ms' # we need a second copy for the hanning smoothed cases 00127 hanningsmooth(vis=clean_inputvis_local_copy2) 00128 00129 00130 # loop over all possible output reference frames 00131 00132 # these are all possible frames: 00133 #frames_to_do = ['TOPO','LSRK', 'LSRD', 'BARY', 'GALACTO', 'LGROUP', 'CMB'] 00134 00135 # the most critical one is CMB (largest freq shift) 00136 #frames_to_do = ['CMB'] 00137 00138 # in order to shorten the test, leave out LSRD, GALACTO, and TOPO 00139 frames_to_do = ['LGROUP', 'LSRK', 'BARY', 'CMB'] 00140 #frames_to_do = ['LSRK'] 00141 00142 for frame in frames_to_do: 00143 00144 restfrq = 22.235080E9 00145 restfreqstr = str(restfrq)+'Hz' 00146 00147 ### frequency mode 00148 00149 outvis = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_freq.ms' 00150 os.system('rm -rf '+outvis) 00151 00152 casalog.post(outvis, 'INFO') 00153 00154 cvel(vis=dataset_name, outputvis=outvis, 00155 mode='frequency',nchan=50, 00156 start=freqmodestart[frame], 00157 width=freqmodewidth[frame], 00158 interpolation='linear', 00159 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00160 outframe=frame, 00161 hanning = dohanning[frame]) 00162 00163 invis = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_freq.ms' 00164 iname = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_freq_clean' 00165 os.system('rm -rf '+iname+'.*') 00166 00167 casalog.post(iname, 'INFO') 00168 00169 clean(vis=invis, 00170 imagename=iname, 00171 field='',spw='', 00172 cell=[0.01,0.01],imsize=[256,256], 00173 stokes='I', 00174 mode='frequency',nchan=50, 00175 start=freqmodestart[frame], 00176 width=freqmodewidth[frame], 00177 interpolation='linear', 00178 psfmode='clark',imagermode='csclean', 00179 scaletype='SAULT', 00180 niter=0,threshold='1.5mJy', 00181 restfreq=restfreqstr, 00182 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00183 mask='', 00184 weighting='briggs', 00185 interactive=F, 00186 minpb=0.3,pbcor=F) 00187 00188 cvel_imstats['frequency'][peakchan][frame] = imstat(iname+'.image', box=testregion, chans=peakchan) 00189 cvel_imstats['frequency'][otherchan1][frame] = imstat(iname+'.image', box=testregion, chans=otherchan1) 00190 cvel_imstats['frequency'][otherchan2][frame] = imstat(iname+'.image', box=testregion, chans=otherchan2) 00191 cvel_imstats['frequency'][otherchan3][frame] = imstat(iname+'.image', box=testregion, chans=otherchan3) 00192 cvel_imstats['frequency'][otherchan4][frame] = imstat(iname+'.image', box=testregion, chans=otherchan4) 00193 00194 iname = 'ANTEN_ALLRE_sort_hann_'+frame+'_freq_clean' 00195 os.system('rm -rf '+iname+'.*') 00196 00197 casalog.post(iname, 'INFO') 00198 00199 cvis = clean_inputvis_local_copy 00200 if(dohanning[frame]): 00201 casalog.post('Will Hanning smooth before cleaning ...', 'INFO') 00202 cvis = clean_inputvis_local_copy2 00203 00204 clean(vis=cvis, 00205 imagename=iname, 00206 field='', spw='', 00207 cell=[0.01,0.01],imsize=[256,256], 00208 stokes='I', 00209 mode='frequency',nchan=50, 00210 start=freqmodestart[frame], 00211 width=freqmodewidth[frame], 00212 outframe=frame, 00213 interpolation='linear', 00214 psfmode='clark',imagermode='csclean', 00215 scaletype='SAULT', 00216 niter=0,threshold='1.5mJy', 00217 restfreq=restfreqstr, 00218 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00219 mask='', 00220 weighting='briggs', 00221 interactive=F, 00222 minpb=0.3,pbcor=F) 00223 00224 cleanonly_imstats['frequency'][peakchan][frame] = imstat(iname+'.image', box=testregion, chans=peakchan) 00225 cleanonly_imstats['frequency'][otherchan1][frame] = imstat(iname+'.image', box=testregion, chans=otherchan1) 00226 cleanonly_imstats['frequency'][otherchan2][frame] = imstat(iname+'.image', box=testregion, chans=otherchan2) 00227 cleanonly_imstats['frequency'][otherchan3][frame] = imstat(iname+'.image', box=testregion, chans=otherchan3) 00228 cleanonly_imstats['frequency'][otherchan4][frame] = imstat(iname+'.image', box=testregion, chans=otherchan4) 00229 00230 00231 #### velocity mode (radio) 00232 00233 f1 = qa.quantity(freqmodestart[frame])['value'] 00234 f2 = f1+qa.quantity(freqmodewidth[frame])['value'] 00235 00236 vrads = (restfrq-f1)/restfrq * 2.99792E8 00237 vradstart = str(vrads)+'m/s' 00238 vradw = (restfrq-f2)/restfrq * 2.99792E8 - vrads 00239 vradwidth = str(vradw)+'m/s' 00240 00241 outvis = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_vrad.ms' 00242 os.system('rm -rf '+outvis) 00243 00244 casalog.post(outvis, 'INFO') 00245 00246 cvel(vis=dataset_name,outputvis=outvis, 00247 mode='velocity',nchan=49, 00248 start=vradstart, 00249 width=vradwidth, 00250 interpolation='fftshift', 00251 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00252 restfreq=restfreqstr, 00253 outframe=frame, 00254 hanning=dohanning[frame]) 00255 00256 invis = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_vrad.ms' 00257 iname = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_vrad_clean' 00258 os.system('rm -rf '+iname+'.*') 00259 00260 casalog.post(iname, 'INFO') 00261 00262 clean(vis=invis, 00263 imagename=iname, 00264 field='',spw='', 00265 cell=[0.01,0.01],imsize=[256,256], 00266 stokes='I', 00267 mode='velocity',nchan=49, 00268 start=vradstart, 00269 width=vradwidth, 00270 interpolation='linear', 00271 psfmode='clark',imagermode='csclean', 00272 scaletype='SAULT', 00273 niter=0,threshold='1.5mJy', 00274 restfreq=restfreqstr, 00275 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00276 mask='', 00277 weighting='briggs', 00278 interactive=F, 00279 minpb=0.3,pbcor=F) 00280 00281 cvel_imstats['radio velocity'][peakchan][frame] = imstat(iname+'.image', box=testregion, chans=peakchan) 00282 cvel_imstats['radio velocity'][otherchan1][frame] = imstat(iname+'.image', box=testregion, chans=otherchan1) 00283 cvel_imstats['radio velocity'][otherchan2][frame] = imstat(iname+'.image', box=testregion, chans=otherchan2) 00284 cvel_imstats['radio velocity'][otherchan3][frame] = imstat(iname+'.image', box=testregion, chans=otherchan3) 00285 cvel_imstats['radio velocity'][otherchan4][frame] = imstat(iname+'.image', box=testregion, chans=otherchan4) 00286 00287 iname = 'ANTEN_ALLRE_sort_hann_'+frame+'_vrad_clean' 00288 os.system('rm -rf '+iname+'.*') 00289 00290 casalog.post(iname, 'INFO') 00291 00292 cvis = clean_inputvis_local_copy 00293 if(dohanning[frame]): 00294 casalog.post('Will Hanning smooth before cleaning ...', 'INFO') 00295 cvis = clean_inputvis_local_copy2 00296 00297 clean(vis=cvis, 00298 imagename=iname, 00299 field='', spw='', 00300 cell=[0.01,0.01],imsize=[256,256], 00301 stokes='I', 00302 mode='velocity',nchan=49, 00303 start=vradstart, 00304 width=vradwidth, 00305 outframe=frame, 00306 interpolation='linear', 00307 psfmode='clark',imagermode='csclean', 00308 scaletype='SAULT', 00309 niter=0,threshold='1.5mJy', 00310 restfreq=restfreqstr, 00311 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00312 mask='', 00313 weighting='briggs', 00314 interactive=F, 00315 minpb=0.3,pbcor=F) 00316 00317 cleanonly_imstats['radio velocity'][peakchan][frame] = imstat(iname+'.image', box=testregion, chans=peakchan) 00318 cleanonly_imstats['radio velocity'][otherchan1][frame] = imstat(iname+'.image', box=testregion, chans=otherchan1) 00319 cleanonly_imstats['radio velocity'][otherchan2][frame] = imstat(iname+'.image', box=testregion, chans=otherchan2) 00320 cleanonly_imstats['radio velocity'][otherchan3][frame] = imstat(iname+'.image', box=testregion, chans=otherchan3) 00321 cleanonly_imstats['radio velocity'][otherchan4][frame] = imstat(iname+'.image', box=testregion, chans=otherchan4) 00322 00323 #### velocity mode (optical) 00324 00325 lambda0 = 2.99792E8/restfrq 00326 lambda1 = 2.99792E8/f1 00327 lambda2 = 2.99792E8/f2 00328 vopts = (lambda1-lambda0)/lambda0 * 2.99792E8 00329 voptw = (lambda2-lambda0)/lambda0 * 2.99792E8 - vopts 00330 voptstart = str(vopts)+'m/s' 00331 voptwidth = str(voptw)+'m/s' 00332 00333 outvis = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_vopt.ms' 00334 os.system('rm -rf '+outvis) 00335 00336 casalog.post(outvis, 'INFO') 00337 00338 cvel(vis=dataset_name, outputvis=outvis, 00339 mode='velocity',nchan=49, 00340 start=voptstart, 00341 width=voptwidth, 00342 interpolation='linear', 00343 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00344 restfreq=restfreqstr, 00345 outframe=frame, 00346 veltype='optical', 00347 hanning=dohanning[frame]) 00348 00349 invis = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_vopt.ms' 00350 iname = 'ANTEN_ALLRE_sort_hann_'+frame+'_cvel_vopt_clean' 00351 os.system('rm -rf '+iname+'.*') 00352 00353 casalog.post(iname, 'INFO') 00354 00355 clean(vis=invis, 00356 imagename=iname, 00357 field='',spw='', 00358 cell=[0.01,0.01],imsize=[256,256], 00359 stokes='I', 00360 mode='velocity',nchan=49, 00361 start=voptstart, 00362 width=voptwidth, 00363 interpolation='linear', 00364 psfmode='clark',imagermode='csclean', 00365 scaletype='SAULT', 00366 niter=0,threshold='1.5mJy', 00367 restfreq=restfreqstr, 00368 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00369 mask='', 00370 weighting='briggs', 00371 interactive=F, 00372 minpb=0.3,pbcor=F, 00373 veltype='optical') 00374 00375 cvel_imstats['optical velocity'][peakchan][frame] = imstat(iname+'.image', box=testregion, chans=peakchan) 00376 cvel_imstats['optical velocity'][otherchan1][frame] = imstat(iname+'.image', box=testregion, chans=otherchan1) 00377 cvel_imstats['optical velocity'][otherchan2][frame] = imstat(iname+'.image', box=testregion, chans=otherchan2) 00378 cvel_imstats['optical velocity'][otherchan3][frame] = imstat(iname+'.image', box=testregion, chans=otherchan3) 00379 cvel_imstats['optical velocity'][otherchan4][frame] = imstat(iname+'.image', box=testregion, chans=otherchan4) 00380 00381 00382 iname = 'ANTEN_ALLRE_sort_hann_'+frame+'_vopt_clean' 00383 os.system('rm -rf '+iname+'.*') 00384 00385 casalog.post(iname, 'INFO') 00386 00387 cvis = clean_inputvis_local_copy 00388 if(dohanning[frame]): 00389 casalog.post('Will Hanning smooth before cleaning ...', 'INFO') 00390 cvis = clean_inputvis_local_copy2 00391 00392 clean(vis=cvis, 00393 imagename=iname, 00394 field='', spw='', 00395 cell=[0.01,0.01],imsize=[256,256], 00396 stokes='I', 00397 mode='velocity',nchan=49, 00398 start=voptstart, 00399 width=voptwidth, 00400 outframe=frame, 00401 interpolation='linear', 00402 psfmode='clark',imagermode='csclean', 00403 scaletype='SAULT', 00404 niter=0,threshold='1.5mJy', 00405 restfreq=restfreqstr, 00406 phasecenter='J2000 12h01m53.13s -18d53m09.8s', 00407 mask='', 00408 weighting='briggs', 00409 interactive=F, 00410 minpb=0.3,pbcor=F, 00411 veltype='optical') 00412 00413 cleanonly_imstats['optical velocity'][peakchan][frame] = imstat(iname+'.image', box=testregion, chans=peakchan) 00414 cleanonly_imstats['optical velocity'][otherchan1][frame] = imstat(iname+'.image', box=testregion, chans=otherchan1) 00415 cleanonly_imstats['optical velocity'][otherchan2][frame] = imstat(iname+'.image', box=testregion, chans=otherchan2) 00416 cleanonly_imstats['optical velocity'][otherchan3][frame] = imstat(iname+'.image', box=testregion, chans=otherchan3) 00417 cleanonly_imstats['optical velocity'][otherchan4][frame] = imstat(iname+'.image', box=testregion, chans=otherchan4) 00418 00419 # end loop over frames 00420 00421 # Analysis 00422 00423 passed = True 00424 # normal values 00425 #tolerance = 0.0013 00426 #avtolerance = 0.00035 00427 # temporarily relaxed values, should be reduced to the normal ones after March 2011 00428 tolerance = 0.012 00429 avtolerance = 0.0039 00430 numpoints = 0. 00431 avdev = 0. 00432 maxdev = 0. 00433 maxdevat = " " 00434 problems = 0 00435 for frame in frames_to_do: 00436 for mode in ['frequency', 'radio velocity', 'optical velocity']: 00437 for chan in mode_imstats.keys(): 00438 print "Testing ", frame, ", ", mode, ", Hanning ", dohanning[frame], ", box ", testregion, ", channel ", chan, " ..." 00439 00440 isok = true 00441 c1 = 0. 00442 c2 = 0. 00443 s1 = 0. 00444 s2 = 0. 00445 try: 00446 c1 = cleanonly_imstats[mode][chan][frame]['max'] 00447 c2 = cvel_imstats[mode][chan][frame]['max'] 00448 s1 = cleanonly_imstats[mode][chan][frame]['maxposf'] 00449 s2 = cvel_imstats[mode][chan][frame]['maxposf'] 00450 except: 00451 isok = False 00452 passed = False 00453 problems +=1 00454 print "Error: this subtest failed" 00455 00456 if isok: 00457 if(abs(c1-c2) > maxdev): 00458 maxdev = abs(c1-c2) 00459 maxdevat = mode+" mode for output frame "+frame\ 00460 +":\n cvel+clean finds max flux in channel "+str(chan)+" to be "+str(c2)\ 00461 +"\n clean-only finds max flux in channel "+str(chan)+" to be "+str(c1) 00462 if(not isnear(c1,c2, tolerance)): 00463 print " ** Problem in ", mode, " mode for output frame ", frame, ":" 00464 print " cvel+clean finds max flux in channel ", chan, " to be ", c2 00465 print " clean-only finds max flux in channel ", chan, " to be ", c1 00466 passed = False 00467 isok = False 00468 problems +=1 00469 00470 avdev += abs(c1-c2) 00471 numpoints += 1. 00472 00473 if(not s1 == s2): 00474 print " ** Problem in ", mode, " mode for output frame ", frame, ":" 00475 print " cvel+clean finds world coordinates for channel ", chan, " to be ", s2 00476 print " clean-only finds world coordinates for channel ", chan, " to be ", s1 00477 passed = False 00478 isok = False 00479 problems +=1 00480 else: 00481 print " World coordinates identical == ", s2 00482 00483 if isok: 00484 print "... OK" 00485 00486 if(numpoints > 0.): 00487 avdev = avdev/numpoints 00488 print numpoints, " spectral points compared, average deviation = ", avdev, " Jy" 00489 if(avdev>avtolerance): 00490 passed = False 00491 print " ** Problem: average deviation too large. Expected is value < ", avtolerance 00492 problems += 1 00493 print " maximum deviation = ", maxdev, " in ", maxdevat 00494 00495 if passed: 00496 print '' 00497 print 'Regression PASSED' 00498 print '' 00499 else: 00500 print "Execution successful but found ", problems, " issues in analysis of results." 00501 print '' 00502 print 'Regression FAILED' 00503 print '' 00504 raise
1.8.0
