|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################################################# 00002 # # 00003 # Test Name: # 00004 # Regression/Benchmarking Script for 3c391 polarization imaging # 00005 # # 00006 # The script is based on EVLA continuum tutorial on 3c391 to test # 00007 # polarization imaging. # 00008 # o starting from calibrated MS # 00009 # o currently only QU imaging with clarkstokes mode is performed # 00010 # mainly to test CAS-2219 fix # 00011 # o does multiscale clean # 00012 # Reference data: # 00013 # o if savecurstats is defined and equal to True, image statistics will # 00014 # saved on a pickle file (this can be used to store to the data # 00015 # repository as a reference data in future regression run. # 00016 # o if pickle file with an appropriate name exist in data directory or # 00017 # current working directoy it try to use the file as reference # 00018 # otherwise use values stored in this script. # 00019 # # 00020 ############################################################################# 00021 import os 00022 import time 00023 from numpy import * 00024 import shutil 00025 import regression_utility as tstutl 00026 00027 #------------------------------------- 00028 # reg test function 00029 def difftest(diffval,tol): 00030 if diffval<tol: 00031 passorfail = 'PASSED' 00032 else: 00033 passorfail = 'FAILED' 00034 return passorfail, diffval 00035 00036 #------------------------------------- 00037 00038 # use Pickle? 00039 # savecurstats will store the current 00040 # regression results in pickle file 00041 loc = locals() 00042 if not loc.has_key('savecurstats'): 00043 savecurstats=False 00044 if savecurstats: 00045 import cPickle 00046 if not loc.has_key('savelog'): 00047 savelog=False 00048 if savelog or savecurstats: 00049 import datetime 00050 00051 # data path 00052 pathname=os.environ.get('CASAPATH').split()[0] 00053 rootdatapath = pathname+'/data/regression/3C391/' 00054 #msdata='3c391_ctm_mosaic_spw0_selfcal.ms' 00055 # Use shorten data: 1 full mosaic sequence (1 sequence each field) 00056 msdata='3c391_ctm_mosaic_spw0_seq1.ms' 00057 00058 # mask (region file) 00059 msk='3c391qu.rgn' 00060 # 00061 if os.path.exists(rootdatapath+msdata): 00062 datapath=rootdatapath+msdata 00063 else: 00064 datapath='./'+msdata 00065 if os.path.exists(rootdatapath+msk): 00066 maskpath=rootdatapath+msk 00067 else: 00068 maskpath='./'+msk 00069 00070 00071 #test working directory 00072 testdir = '3c391polim_regression' 00073 00074 # clean up previous run 00075 tstutl.cleanup(testdir) 00076 00077 # copying the data to testdir 00078 tstutl.maketestdir(testdir) 00079 print "Test working directory: ", testdir 00080 try: 00081 print "copying data from "+datapath 00082 shutil.copytree(datapath,testdir+'/'+msdata) 00083 print "copying mask image from "+maskpath 00084 shutil.copyfile(maskpath,testdir+'/'+msk) 00085 except Exception,e: 00086 tstutl.note("Failed with:"+str(e)) 00087 raise e 00088 00089 ################## 00090 # reference data # 00091 ################## 00092 refstats={} 00093 # 00094 # Q image 00095 #*.image 00096 # on rishi (casapy-test 3.1.0 r13281) 00097 refstats['qim'] = {'blc': array([0, 0, 0, 0], dtype=int32), 00098 'blcf': '18:50:12.253, -01.07.40.561, Q, 4.599e+09Hz', 00099 'flux': array([-0.02223448]), 00100 'max': array([ 0.00106219]), 00101 'maxpos': array([339, 268, 0, 0], dtype=int32), 00102 'maxposf': '18:49:15.743, -00.56.30.579, Q, 4.599e+09Hz', 00103 'mean': array([ -9.41242758e-06]), 00104 'medabsdevmed': array([ 3.70866292e-05]), 00105 'median': array([ -1.42126601e-06]), 00106 'min': array([-0.00229603]), 00107 'minpos': array([334, 300, 0, 0], dtype=int32), 00108 'minposf': '18:49:16.576, -00.55.10.579, Q, 4.599e+09Hz', 00109 'npts': array([ 105531.]), 00110 'quartile': array([ 7.41678959e-05]), 00111 'rms': array([ 0.00012677]), 00112 'sigma': array([ 0.00012642]), 00113 'sum': array([-0.9933029]), 00114 'sumsq': array([ 0.00169583]), 00115 'trc': array([575, 575, 0, 0], dtype=int32), 00116 'trcf': '18:48:36.407, -00.43.43.058, Q, 4.599e+09Hz'} 00117 # .residual 00118 refstats['qres'] = {'blc': array([0, 0, 0, 0], dtype=int32), 00119 'blcf': '18:50:12.253, -01.07.40.561, Q, 4.599e+09Hz', 00120 'max': array([ 0.00100225]), 00121 'maxpos': array([250, 92, 0, 0], dtype=int32), 00122 'maxposf': '18:49:30.578, -01.03.50.580, Q, 4.599e+09Hz', 00123 'mean': array([ 3.80984632e-09]), 00124 'medabsdevmed': array([ 5.04056925e-06]), 00125 'median': array([ 0.]), 00126 'min': array([-0.0007955]), 00127 'minpos': array([490, 306, 0, 0], dtype=int32), 00128 'minposf': '18:48:50.573, -00.54.55.570, Q, 4.599e+09Hz', 00129 'npts': array([ 331776.]), 00130 'quartile': array([ 1.00784273e-05]), 00131 'rms': array([ 3.89150882e-05]), 00132 'sigma': array([ 3.89151474e-05]), 00133 'sum': array([ 0.00126402]), 00134 'sumsq': array([ 0.00050244]), 00135 'trc': array([575, 575, 0, 0], dtype=int32), 00136 'trcf': '18:48:36.407, -00.43.43.058, Q, 4.599e+09Hz'} 00137 # U image stats on rishi 00138 #*.image 00139 refstats['uim'] = {'blc': array([0, 0, 1, 0], dtype=int32), 00140 'blcf': '18:50:12.253, -01.07.40.561, U, 4.599e+09Hz', 00141 'flux': array([-0.00019476]), 00142 'max': array([ 0.00213059]), 00143 'maxpos': array([283, 336, 1, 0], dtype=int32), 00144 'maxposf': '18:49:25.077, -00.53.40.580, U, 4.599e+09Hz', 00145 'mean': array([ -8.24477105e-08]), 00146 'medabsdevmed': array([ 4.60967676e-05]), 00147 'median': array([ -8.16736417e-07]), 00148 'min': array([-0.00212354]), 00149 'minpos': array([330, 302, 1, 0], dtype=int32), 00150 'minposf': '18:49:17.243, -00.55.05.580, U, 4.599e+09Hz', 00151 'npts': array([ 105531.]), 00152 'quartile': array([ 9.21983010e-05]), 00153 'rms': array([ 0.00013302]), 00154 'sigma': array([ 0.00013302]), 00155 'sum': array([-0.00870079]), 00156 'sumsq': array([ 0.00186722]), 00157 'trc': array([575, 575, 1, 0], dtype=int32), 00158 'trcf': '18:48:36.407, -00.43.43.058, U, 4.599e+09Hz'} 00159 00160 #*.residual 00161 refstats['ures'] = {'blc': array([0, 0, 1, 0], dtype=int32), 00162 'blcf': '18:50:12.253, -01.07.40.561, U, 4.599e+09Hz', 00163 'max': array([ 0.0009961]), 00164 'maxpos': array([355, 92, 1, 0], dtype=int32), 00165 'maxposf': '18:49:13.075, -01.03.50.579, U, 4.599e+09Hz', 00166 'mean': array([ -1.25753895e-07]), 00167 'medabsdevmed': array([ 5.80839878e-06]), 00168 'median': array([ 0.]), 00169 'min': array([-0.0008242]), 00170 'minpos': array([270, 484, 1, 0], dtype=int32), 00171 'minposf': '18:49:27.244, -00.47.30.579, U, 4.599e+09Hz', 00172 'npts': array([ 331776.]), 00173 'quartile': array([ 1.16159445e-05]), 00174 'rms': array([ 4.75954876e-05]), 00175 'sigma': array([ 4.75953917e-05]), 00176 'sum': array([-0.04172212]), 00177 'sumsq': array([ 0.00075158]), 00178 'trc': array([575, 575, 1, 0], dtype=int32), 00179 'trcf': '18:48:36.407, -00.43.43.058, U, 4.599e+09Hz'} 00180 00181 # regression state intialization 00182 regstate = True 00183 00184 ############# 00185 #regression # 00186 ############# 00187 try: 00188 startTime = time.time() 00189 startProc = time.clock() 00190 00191 # clean 00192 print "***Imaging***" 00193 quimagename='3c391_ctm_spw0_QU_selfcal' 00194 00195 quimagepath=testdir+'/'+quimagename 00196 ret=clean(vis=testdir+'/'+msdata,imagename=quimagepath, 00197 field='',spw='', 00198 mode='mfs', 00199 #niter=5000, 00200 niter=8000, 00201 gain=0.1,threshold='0.0mJy', 00202 psfmode='clarkstokes', 00203 imagermode='mosaic',ftmachine='mosaic', 00204 multiscale=[0, 6, 18, 54],smallscalebias=0.9, 00205 #interactive=True, 00206 mask=testdir+'/'+msk, 00207 imsize=[576,576],cell=['2.5arcsec','2.5arcsec'], 00208 stokes='QU', 00209 weighting='briggs',robust=0.0, 00210 usescratch=False) 00211 00212 endTime = time.time() 00213 endProc = time.clock() 00214 except: 00215 print "Clean execution failed" 00216 00217 else: 00218 # analysis 00219 print "***Analyze the results***" 00220 curstats={} 00221 imagename=quimagepath+'.image' 00222 curstats['qim']=imstat(imagename=imagename, stokes='Q') 00223 curstats['uim']=imstat(imagename=imagename, stokes='U') 00224 imagename=quimagepath+'.residual' 00225 curstats['qres']=imstat(imagename=imagename, stokes='Q') 00226 curstats['ures']=imstat(imagename=imagename, stokes='U') 00227 00228 savelog=False 00229 if savelog: 00230 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00231 outfile='3c391polimg.'+datestring+'.log' 00232 logfile=open(outfile,'w') 00233 else: 00234 logfile=sys.stdout 00235 00236 # get reference data if pickle file is exist 00237 refdataused='reference values defined inside the script' 00238 refpath=rootdatapath+'reference/3C391linpol_regrun.pickle' 00239 if os.path.isfile(refpath): 00240 f = open(refpath) 00241 import cPickle 00242 try: 00243 refarch=cPickle.load(f) 00244 refstattypes=cPickle.load(f) 00245 for sttype in refstattypes: 00246 refstats[sttype]=cPickle.load(f) 00247 refdataused=refpath 00248 except (EOFError, cPickle.UnpicklingError): 00249 raise e 00250 f.close() 00251 00252 # tolerance 5% 00253 tol = 0.05 00254 # as of Oct 21, 2010: it does not seem to pass with 5% for U residual test 00255 # on active (its ok for prerelease) on 32bit linux. So relax the value a bit 00256 # for now. 00257 tolres = 0.085 00258 00259 # list of deviation evaluation 00260 # currently tests max flux density in cleaned image and sigma in residual image 00261 diffqmaxfrac=abs((curstats['qim']['max'][0] - refstats['qim']['max'][0])/refstats['qim']['max'][0]) 00262 diffumaxfrac=abs((curstats['uim']['max'][0] - refstats['uim']['max'][0])/refstats['uim']['max'][0]) 00263 diffqsigfrac=abs((curstats['qres']['sigma'][0] - refstats['qres']['sigma'][0])/refstats['qres']['sigma'][0]) 00264 diffusigfrac=abs((curstats['ures']['sigma'][0] - refstats['ures']['sigma'][0])/refstats['ures']['sigma'][0]) 00265 00266 # acutual tests 00267 tests={} 00268 tests['qmaxdiff']=difftest(diffqmaxfrac,tol) 00269 tests['qsigdiff']=difftest(diffqsigfrac,tolres) 00270 tests['umaxdiff']=difftest(diffumaxfrac,tol) 00271 tests['usigdiff']=difftest(diffusigfrac,tolres) 00272 00273 # results logging 00274 print "***Logging the reulsts***" 00275 failcnts=0 00276 for k,v in tests.items(): 00277 if v=='FAILED': 00278 regstate=False 00279 failcnts+=1 00280 if regstate: 00281 msg='PASSED' 00282 else: 00283 msg='FAILED (ntest failed=%s)' % failcnts 00284 00285 print >>logfile,'' 00286 print >>logfile,'' 00287 print >>logfile,'********************************* Data Summary *********************************' 00288 print >>logfile,'* ' 00289 print >>logfile,'Observation: EVLA' 00290 print >>logfile,'Data records: 54582 Total integration time = 2085.5 seconds' 00291 print >>logfile,'Observed from 24-Apr-2010/08:24:53.0 to 24-Apr-2010/08:59:38.5 (UTC)' 00292 print >>logfile,'ObservationID = 0 ArrayID = 0' 00293 print >>logfile,'Date Timerange (UTC) Scan FldId FieldName nVis Int(s) SpwIds' 00294 print >>logfile,'24-Apr-2010/08:24:53.0 - 08:29:43.0 5 0 3C391 C1 7590 9.83 [0]' 00295 print >>logfile,' 08:29:43.0 - 08:34:43.0 6 1 3C391 C2 7821 9.65 [0]' 00296 print >>logfile,' 08:34:43.0 - 08:39:43.0 7 2 3C391 C3 7821 9.66 [0]' 00297 print >>logfile,' 08:39:43.0 - 08:44:43.0 8 3 3C391 C4 7821 9.61 [0]' 00298 print >>logfile,' 08:44:43.0 - 08:49:43.0 9 4 3C391 C5 7843 9.62 [0]' 00299 print >>logfile,' 08:49:43.0 - 08:54:43.0 10 5 3C391 C6 7843 9.58 [0]' 00300 print >>logfile,' 08:54:43.0 - 08:59:38.5 11 6 3C391 C7 7843 9.58 [0]' 00301 print >>logfile,' (nVis = Total number of time/baseline visibilities per scan)' 00302 print >>logfile,'Fields: 7' 00303 print >>logfile,' ID Code Name RA Decl Epoch SrcId nVis' 00304 print >>logfile,' 0 NONE 3C391 C1 18:49:24.2440 -00.55.40.5800 J2000 0 7590' 00305 print >>logfile,' 1 NONE 3C391 C2 18:49:29.1490 -00.57.48.0000 J2000 1 7821' 00306 print >>logfile,' 2 NONE 3C391 C3 18:49:19.3390 -00.57.48.0000 J2000 2 7821' 00307 print >>logfile,' 3 NONE 3C391 C4 18:49:14.4340 -00.55.40.5800 J2000 3 7821' 00308 print >>logfile,' 4 NONE 3C391 C5 18:49:19.3390 -00.53.33.1600 J2000 4 7843' 00309 print >>logfile,' 5 NONE 3C391 C6 18:49:29.1490 -00.53.33.1600 J2000 5 7843' 00310 print >>logfile,' 6 NONE 3C391 C7 18:49:34.0540 -00.55.40.5800 J2000 6 7843' 00311 print >>logfile,' (nVis = Total number of time/baseline visibilities per field)' 00312 print >>logfile,'Spectral Windows: (1 unique spectral windows and 1 unique polarization setups)' 00313 print >>logfile,' SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs' 00314 print >>logfile,' 0 64 TOPO 4536 2000 128000 4536 RR RL LR LL ' 00315 print >>logfile,'********************************************************************************' 00316 print >>logfile,' ' 00317 print >>logfile,'******************************** Regression ************************************' 00318 print >>logfile,'* *' 00319 print >>logfile,'Q image test' 00320 print >>logfile,'Peak flux density: %-10.5g (expected: %-10.5g, frac_diff:%-10.5g) %s' \ 00321 % (curstats['qim']['max'][0], refstats['qim']['max'][0],\ 00322 tests['qmaxdiff'][1], tests['qmaxdiff'][0]) 00323 print >>logfile,'(tolerance=%s)' % tol 00324 print >>logfile,'Residual sigma: %-10.5g (expected: %-10.5g, frac_diff:%-10.5g) %s ' \ 00325 % (curstats['qres']['sigma'][0], refstats['qres']['sigma'][0],\ 00326 tests['qsigdiff'][1],tests['qsigdiff'][0]) 00327 print >>logfile,'(tolerance=%s)' % tolres 00328 print >>logfile,'------------------------------------------------------------' 00329 print >>logfile,'U image test ' 00330 print >>logfile,'Peak flux density: %-10.5g (expected: %-10.5g, frac_diff:%-10.5g) %s' \ 00331 % (curstats['uim']['max'][0], refstats['uim']['max'][0], \ 00332 tests['umaxdiff'][1], tests['umaxdiff'][0]) 00333 print >>logfile,'(tolerance=%s)' % tol 00334 print >>logfile,'Residual sigma: %-10.5g (expected: %-10.5g, frac_diff:%-10.5g) %s ' \ 00335 % (curstats['ures']['sigma'][0], refstats['ures']['sigma'][0],\ 00336 tests['usigdiff'][1],tests['usigdiff'][0]) 00337 print >>logfile,'(tolerance=%s)' % tolres 00338 print >>logfile,' ' 00339 print >>logfile,' ** reference data used ** ' 00340 print >>logfile, refdataused 00341 print >>logfile,' ============================================================================== ' 00342 print >>logfile,'Regression %s ' % msg 00343 print >>logfile,'* *' 00344 print >>logfile,'********************************************************************************' 00345 print >>logfile,' ' 00346 print >>logfile,'******************** Benchmarking **************************' 00347 print >>logfile,'* *' 00348 print >>logfile,'Total wall clock time was %13.3f: ' % (endTime - startTime) 00349 print >>logfile,'Total CPU time was %13.3f: ' % (endProc - startProc) 00350 print >>logfile,'Processing rate MB/s was %13.3f: ' % (329./(endTime - startTime)) 00351 print >>logfile,'* *' 00352 print >>logfile,'************************************************************' 00353 if savelog: logfile.close() 00354 00355 # save the regression results (statistics) in pickle 00356 # for reference data 00357 if savecurstats: 00358 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00359 pickfile = '3c391polim.regression.'+datestring+'.pickle' 00360 f = open(pickfile,'w') 00361 casapath=os.environ['CASAPATH'].split() 00362 if os.environ.has_key('HOSTNAME'): 00363 hostname=os.environ['HOSTNAME'] 00364 elif len(casapath) >3: 00365 hostname=casapath[3] 00366 else: 00367 hostname='unknown' 00368 # check the environment that the script is running 00369 if casapath[1] =='darwin': 00370 # darwin 00371 # use sw_vers to get OSX version ? 00372 os.system('/usr/bin/sw_vers >/tmp/osxverinfo.txt') 00373 fver=open('/tmp/osxverinfo.txt') 00374 fver.readline() 00375 verinfo=fver.readline().split(':')[1].strip() 00376 arch="running on %s (Mac OS X %s)" % (hostname, verinfo) 00377 elif pathname.find('lib64')>-1: 00378 # 64bit linux 00379 arch="running on %s (64bit linux)" % hostname 00380 else: 00381 # assume 32bit linux 00382 arch ="running on %s (32bit linux)" % hostname 00383 #Pickling 00384 cPickle.dump(arch,f) 00385 # stats are stored in this order 00386 statslist=['qim','qres','uim','ures'] 00387 cPickle.dump(statslist,f) 00388 for stats in statslist: 00389 cPickle.dump(curstats[stats],f) 00390 f.close() 00391 00392 # just to be informative 00393 if savelog: 00394 print "Regression %s: " % msg 00395 print "log output: %s" % outfile
1.8.0
