|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ################################################################################ 00002 # # 00003 # Regression/Benchmarking Script for wideband imaging # 00004 # # 00005 ################################################################################ 00006 # # 00007 # (1) MS-MFS with nterms=3 on VLA_multifrequency_3C286.ms # 00008 # - checks output residual rms and total power # 00009 # - checks output peak intensity, spectral index and spectral curvature # 00010 # 00011 # (2) Post-deconvolution PB-correction, on the output of (1) 00012 # - checks output peak corrected intensity, spectral index and curvature. 00013 # # 00014 ################################################################################ 00015 # # 00016 # More tests that will appear here in the future : # 00017 # # 00018 # (2) MS-MFS with wide-band primary-beam correction # 00019 # (3) MS-MFS on extended emission : DONE # 00020 # (4) MS-MFS with mosaicing # 00021 # # 00022 ################################################################################ 00023 00024 import time 00025 import os 00026 00027 # Data : VLA_multifrequency_3C286.ms 00028 pathname=os.environ.get('CASAPATH').split()[0] 00029 00030 # Initialize status flag 00031 regstate = True; 00032 00033 # Start timers 00034 startTime=time.time() 00035 startProc=time.clock() 00036 00037 # Mark time 00038 copyTime=time.time() 00039 00040 if(regstate): 00041 # Test (1) : Run the clean task 00042 print '--Image with MS-MFS--' 00043 default('clean') 00044 npix=1024; 00045 ret = clean(vis='VLA_multifrequency_3C286.ms',imagename='reg_3C286',imagermode='', 00046 nterms=3,reffreq='1.4GHz', 00047 niter=50,gain=0.8,threshold='7.0mJy',imsize=[npix,npix], 00048 cell=['2.5arcsec','2.5arcsec'],weighting='briggs',usescratch=False); 00049 00050 # Test (2) : Post-deconvolution PB-correction 00051 print '--Post deconvolution wideband PB correction--' 00052 ret = widebandpbcor(vis='VLA_multifrequency_3C286.ms', imagename='reg_3C286', 00053 nterms=3, threshold='5mJy', action='pbcor', reffreq='1.4GHz', 00054 pbmin=0.2, field='',spwlist=[0, 1, 2, 3, 4, 5, 6], chanlist=[8, 8, 8, 8, 8, 8, 8], 00055 weightlist=[1, 1, 1, 1, 1, 1, 1]) 00056 00057 # Stop timers 00058 endProc=time.clock() 00059 endTime=time.time() 00060 00061 # Start printing info. 00062 import datetime 00063 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00064 outfile='reg_3C286.'+datestring+'.log' 00065 logfile=open(outfile,'w') 00066 00067 # Data summary 00068 print >>logfile,'**************** Regression-test for wide-band imaging *************************' 00069 print >>logfile,'** **' 00070 print >>logfile,'****************** 3C286 wideband data (L-Band) *******************************' 00071 print >>logfile,'' 00072 print >>logfile,'Observation: VLA' 00073 print >>logfile,'Data records: 145600 Total integration time = 35435 seconds' 00074 print >>logfile,' Observed from 28-Apr-2008/00:53:22.5 to 28-Apr-2008/10:43:57.5 (UTC)' 00075 print >>logfile,'Fields: 1' 00076 print >>logfile,' ID Code Name RA Decl Epoch SrcId nVis ' 00077 print >>logfile,' 0 A 1331+305 13:31:08.2879 +30.30.32.9580 J2000 0 145600 ' 00078 print >>logfile,' (nVis = Total number of time/baseline visibilities per field) ' 00079 print >>logfile,'Spectral Windows: (7 unique spectral windows and 1 unique polarization setups)' 00080 print >>logfile,' SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs ' 00081 print >>logfile,' 0 15 TOPO 1184.0625 1562.5 23437.5 1195 RR LL ' 00082 print >>logfile,' 1 15 TOPO 1301.0625 1562.5 23437.5 1312 RR LL ' 00083 print >>logfile,' 2 15 TOPO 1401.0625 1562.5 23437.5 1412 RR LL ' 00084 print >>logfile,' 3 15 TOPO 1494.0625 1562.5 23437.5 1505 RR LL ' 00085 print >>logfile,' 4 15 TOPO 1676.0625 1562.5 23437.5 1687 RR LL ' 00086 print >>logfile,' 5 15 TOPO 1751.0625 1562.5 23437.5 1762 RR LL ' 00087 print >>logfile,' 6 15 TOPO 1864.0625 1562.5 23437.5 1875 RR LL ' 00088 00089 00090 # Perform the checks 00091 print >>logfile,'' 00092 print >>logfile,'*********************** Comparison of results **********************************' 00093 print >>logfile,'** **' 00094 print >>logfile,'** (1) MS-MFS on the 3C286 field with nterms=3 and reffreq=1.4GHz **' 00095 print >>logfile,'** **' 00096 print >>logfile,'********************************************************************************' 00097 00098 if(not regstate): 00099 print >>logfile,'* Data file VLA_multifrequency_3C286.ms cannot be found'; 00100 else: 00101 # (6) : (V2.5) This is the truth (for active, 19Jan2012) - wrote coefficient residuals to the output residual image, instead of using them only for alpha,beta calculations 00102 # Changes from previous numbers 'active r17725' are mainly 'noise' levels. 00103 correct_sigma = 0.00126019103471; 00104 correct_sumsq = 1.6652231052; 00105 correct_intensity = 14.8404045105; 00106 correct_alpha = -0.471577763557; 00107 correct_beta = -0.124552100897; 00108 ## Added on 13 Sep 2012 00109 correct_pbcor_intensity = 0.234191760421 00110 correct_pbcor_alpha = -0.910139858723 00111 00112 # (5) : (V2.3) This is the truth (for active, 02Mar2010) - included coefficient residuals in alpha/beta calcs. 00113 # Changes from previous numbers 'active r14198' are within the noise, and only for alpha/beta. 00114 #correct_sigma = 0.00095682840; 00115 #correct_sumsq = 0.959992192; 00116 #correct_intensity = 14.84169483; 00117 #correct_alpha = -0.4716707468; 00118 #correct_beta = -0.124309256; 00119 00120 # (4) : (V2.2) This is the truth (for active, 21Jan2010) - removed extra vecpsf0 conv. 00121 # Changes from previous numbers 'active r13845' are within the noise. 00122 #correct_sigma = 0.00095682840; 00123 #correct_sumsq = 0.959992192; 00124 #correct_intensity = 14.84169483; 00125 #correct_alpha = -0.471375882; 00126 #correct_beta = -0.127162337; 00127 00128 # (3) : This is the truth (for active, 20Dec2010) - change from MTLC to MTMC (+ algorithm fiddling). 00129 # Changes from previous numbers for 'active r13787' are within the noise. 00130 #correct_sigma = 0.00095900; 00131 #correct_sumsq = 0.9644402; 00132 #correct_intensity = 14.8406848; 00133 #correct_alpha = -0.47158026; 00134 #correct_beta = -0.12506663; 00135 00136 # (2) : This is the truth (for active, 21Oct2010) - with SB's gridding fixes 00137 #correct_sigma = 0.00099339; 00138 #correct_sumsq = 1.03476342; 00139 #correct_intensity = 14.8406724; 00140 #correct_alpha = -0.4706874; 00141 #correct_beta = -0.12786445; 00142 00143 # (1) : This is the truth (for prerelease, 21Oct2010) - without SB's gridding fixes. 00144 #correct_sigma = 0.0010294; 00145 #correct_sumsq = 1.11118678; 00146 #correct_intensity = 14.838494; 00147 #correct_alpha = -0.47109225; 00148 #correct_beta = -0.12466369; 00149 00150 # Residual rms noise and sum-sq (total power) 00151 if(os.path.exists('reg_3C286.residual.tt0')): 00152 ia.open('reg_3C286.residual.tt0'); 00153 stats = ia.statistics(); 00154 ia.close(); 00155 diff_sigma = abs( (stats['sigma'][0]) - correct_sigma )/correct_sigma; 00156 diff_sumsq = abs( (stats['sumsq'][0]) - correct_sumsq )/correct_sumsq; 00157 if(diff_sigma<0.05): 00158 print >>logfile,'* Passed residual sigma test '; 00159 else: 00160 print >>logfile,'* FAILED residual sigma test ' 00161 regstate = False; 00162 print >>logfile,'-- residual sigma : ' + str((stats['sigma'][0])) + ' (' + str(correct_sigma) + ')'; 00163 if(diff_sumsq<0.05): 00164 print >>logfile,'* Passed residual total-power test '; 00165 else: 00166 print >>logfile,'* FAILED residual total-power test ' 00167 regstate = False 00168 print >>logfile,'-- residual sumsq : ' + str((stats['sumsq'][0])) + ' (' + str(correct_sumsq) + ')'; 00169 else: 00170 print >>logfile,' FAILED : No residual image generated.' 00171 regstate = False; 00172 00173 # Intensity 00174 if(os.path.exists('reg_3C286.image.tt0')): 00175 ia.open('reg_3C286.image.tt0'); 00176 midpix = ia.pixelvalue([npix/2,npix/2]) 00177 ia.close(); 00178 diff_intensity = abs( midpix['value']['value'] - correct_intensity )/ abs(correct_intensity); 00179 if(diff_intensity<0.02): 00180 print >>logfile,'* Passed peak intensity test '; 00181 else: 00182 print >>logfile,'* FAILED peak intensity test ' 00183 regstate = False; 00184 print >>logfile,'-- peak intensity : ' + str(midpix['value']['value']) + ' (' + str(correct_intensity) + ')'; 00185 else: 00186 print >>logfile,'-- FAILED : No intensity map generated'; 00187 regstate = False; 00188 00189 # Alpha 00190 if(os.path.exists('reg_3C286.image.alpha')): 00191 ia.open('reg_3C286.image.alpha'); 00192 midpix = ia.pixelvalue([npix/2,npix/2]) 00193 ia.close(); 00194 diff_alpha = abs( midpix['value']['value'] - correct_alpha )/ abs(correct_alpha); 00195 if(diff_alpha<0.02): 00196 print >>logfile,'* Passed spectral index test '; 00197 else: 00198 print >>logfile,'* FAILED spectral index test ' 00199 regstate = False; 00200 print >>logfile,'-- spectral index : ' + str(midpix['value']['value']) + ' (' + str(correct_alpha) + ')'; 00201 else: 00202 print >>logfile,'-- FAILED : No spectral index map generated'; 00203 regstate = False; 00204 00205 # Beta 00206 if(os.path.exists('reg_3C286.image.beta')): 00207 ia.open('reg_3C286.image.beta'); 00208 midpix = ia.pixelvalue([npix/2,npix/2]) 00209 ia.close(); 00210 diff_beta = abs( midpix['value']['value'] - correct_beta )/ abs(correct_beta); 00211 if(diff_beta<0.02): 00212 print >>logfile,'* Passed spectral curvature test '; 00213 else: 00214 print >>logfile,'* FAILED spectral curvature test ' 00215 regstate = False; 00216 print >>logfile,'-- spectral curvature : ' + str(midpix['value']['value']) + ' (' + str(correct_beta) + ')'; 00217 else: 00218 print >>logfile,'-- FAILED : No spectral curvature map generated'; 00219 regstate = False; 00220 00221 # PB-corrected intensity) 00222 if(os.path.exists('reg_3C286.pbcor.image.tt0')): 00223 ia.open('reg_3C286.pbcor.image.tt0'); 00224 offpix = ia.pixelvalue([304,542]) 00225 ia.close(); 00226 diff_int = abs( offpix['value']['value'] - correct_pbcor_intensity )/ abs(correct_pbcor_intensity); 00227 if(diff_int<0.02): 00228 print >>logfile,'* Passed widebandpbcor intensity test '; 00229 else: 00230 print >>logfile,'* FAILED widebandpbcor intensity test ' 00231 regstate = False; 00232 print >>logfile,'-- pb-corrected intensity : ' + str(offpix['value']['value']) + ' (' + str(correct_pbcor_intensity) + ')'; 00233 else: 00234 print >>logfile,'-- FAILED : No pb-corrected intensity map generated'; 00235 regstate = False; 00236 00237 # PB-corrected alpha (not checking beta) 00238 if(os.path.exists('reg_3C286.pbcor.image.alpha')): 00239 ia.open('reg_3C286.pbcor.image.alpha'); 00240 offpix = ia.pixelvalue([304,542]) 00241 ia.close(); 00242 diff_alpha = abs( offpix['value']['value'] - correct_pbcor_alpha )/ abs(correct_pbcor_alpha); 00243 if(diff_alpha<0.02): 00244 print >>logfile,'* Passed widebandpbcor alpha test '; 00245 else: 00246 print >>logfile,'* FAILED widebandpbcor alpha test ' 00247 regstate = False; 00248 print >>logfile,'-- pb-corrected spectral index : ' + str(offpix['value']['value']) + ' (' + str(correct_pbcor_alpha) + ')'; 00249 else: 00250 print >>logfile,'-- FAILED : No pb-corrected spectral index map generated'; 00251 regstate = False; 00252 00253 00254 # Final verdict 00255 if(regstate): 00256 print >>logfile,'PASSED regression test for wideband-imaging.' 00257 print '' 00258 print 'Regression PASSED' 00259 print '' 00260 else: 00261 print >>logfile,'FAILED regression test for wideband-imaging.' 00262 print '' 00263 print 'Regression FAILED' 00264 print '' 00265 00266 print >>logfile,'' 00267 00268 # Print timing info 00269 print >>logfile,'********************************************************************************' 00270 print >>logfile,'** Benchmarking **' 00271 print >>logfile,'********************************************************************************' 00272 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00273 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00274 print >>logfile,'Processing rate MB/s was: '+str(278./(endTime - startTime)) 00275 print >>logfile,'* Breakdown: *' 00276 print >>logfile,'* copy time was: '+str(copyTime-startTime) 00277 print >>logfile,'* imaging time was: '+str(endTime-copyTime) 00278 print >>logfile,'* *' 00279 print >>logfile,'********************************************************************************' 00280 00281 logfile.close() 00282 00283
1.8.0
