|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ####################################################### 00002 # # 00003 # Regression/Test Script for Stokes-Selection in Imager and Clean # 00004 # # 00005 ####################################################### 00006 # 00007 # (1) Given a 'true' [I,Q,U,V] vector, 00008 # simulate MSs for [RR,RL,LR,LL], [RR,LL], [RR],[LL], [XX,XY,YX,YY], [XX,YY], [XX], [YY]. 00009 # (2) Test all possible options, some valid, some not : ['I','Q','U','V', 'IV', 'QU', 'IQ', 'UV', 00010 # 'IQU', 'IUV', 'IQUV', 'RR', 'LL', 'RL','LR', 'RRLL','RLLR','XX','YY','XY','YX','XXYY','XYYX'] 00011 # (3) For runs that produce an output image, the resulting output (peak pixel value) 00012 # is compared with the 'true' [I,Q,U,V,RR,LL,XX,YY]. 00013 # (4) Counts of 'pass' and 'fail' are, in the end, compared and used to detect regression. 00014 # Note : This script tests both valid and invalid options. 00015 # 00016 # 00017 ###################################################### 00018 # 00019 # More tests that will appear here in the future 00020 # (See JIRA CAS 2587 : https://bugs.nrao.edu/browse/CAS-2587 ) : 00021 # 00022 # (1) Tests with partially flagged data 00023 # (2) Tests with different minor-cycle algorithms 00024 # 00025 ###################################################### 00026 00027 import time 00028 import os 00029 00030 # Init printing info. 00031 import datetime 00032 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00033 outfile='stokesimagingtest.'+datestring+'.log' 00034 logfile=open(outfile,'w') 00035 00036 ##################################################### 00037 #### BEGIN : Helper function definition (tests are run after 'END') #### 00038 ##################################################### 00039 00040 ################################ 00041 ### Simulate a small MS 00042 ################################ 00043 def makeMS(choice='RL_all',stokesvals=[1.0,1.0,0.0,0.0],parentpath='.'): 00044 dirname = parentpath+"/Point" 00045 choicelist = ['RL_all','RL_par','RR_only','LL_only','XY_all','XY_par','XX_only','YY_only']; 00046 if(not (choice in choicelist)): 00047 print >>logfile, 'Cannot find specs for ', choice; 00048 return; 00049 if(choice=='RL_all'): 00050 basename = dirname+"/point_linRL" 00051 msstokes='RR RL LR LL'; 00052 feedtype='perfect R L'; 00053 if(choice=='XY_all'): 00054 basename = dirname+"/point_linXY" 00055 msstokes='XX XY YX YY'; 00056 feedtype='perfect X Y'; 00057 if(choice=='RL_par'): 00058 basename = dirname+"/point_parRL" 00059 msstokes='RR LL'; 00060 feedtype='perfect R L'; 00061 if(choice=='XY_par'): 00062 basename = dirname+"/point_parXY" 00063 msstokes='XX YY'; 00064 feedtype='perfect X Y'; 00065 if(choice=='XX_only'): 00066 basename = dirname+"/point_onlyXX" 00067 msstokes='XX'; 00068 feedtype='perfect X Y'; 00069 if(choice=='YY_only'): 00070 basename = dirname+"/point_onlyYY" 00071 msstokes='YY'; 00072 feedtype='perfect X Y'; 00073 if(choice=='RR_only'): 00074 basename = dirname+"/point_onlyRR" 00075 msstokes='RR'; 00076 feedtype='perfect R L'; 00077 if(choice=='LL_only'): 00078 basename = dirname+"/point_onlyLL" 00079 msstokes='LL'; 00080 feedtype='perfect R L'; 00081 00082 msname = basename + '.ms'; 00083 clname = basename+'.cl'; 00084 00085 if(not os.path.exists(dirname)): 00086 cmd = 'mkdir ' + dirname; 00087 os.system(cmd); 00088 00089 vx = [41.1100006, 134.110001, 268.309998, 439.410004, 644.210022] 00090 vy = [3.51999998, -39.8300018, -102.480003, -182.149994, -277.589996] 00091 vz = [0.25, -0.439999998, -1.46000004, -3.77999997, -5.9000001] 00092 d = [25.0, 25.0, 25.0, 25.0, 25.0] 00093 an = ['VLA1','VLA2','VLA3','VLA4','VLA5']; 00094 nn = len(vx); 00095 x = (vx - (sum(pl.array(vx))/(nn))); 00096 y = (vy - (sum(pl.array(vy))/(nn))); 00097 z = (vz - (sum(pl.array(vz))/(nn))); 00098 00099 sm.open(ms=msname); 00100 sm.setconfig(telescopename='VLA',x=x.tolist(),y=y.tolist(),z=z.tolist(),dishdiameter=d, 00101 mount=['alt-az'], antname=an, 00102 coordsystem='local',referencelocation=me.observatory('VLA')); 00103 sm.setspwindow(spwname="LBand",freq="1.0GHz",deltafreq='500MHz', 00104 freqresolution='2MHz',nchannels=3,stokes=msstokes); 00105 sm.setfeed(mode=feedtype,pol=['']); 00106 ra0="19:59:28.500"; 00107 dec0="+40.44.01.50"; 00108 sm.setfield( sourcename="fake",sourcedirection=me.direction(rf='J2000',v0=ra0,v1=dec0) ); 00109 sm.setlimits(shadowlimit=0.01, elevationlimit='10deg'); 00110 sm.setauto(autocorrwt=0.0); 00111 sm.settimes(integrationtime='2000s', usehourangle=True, 00112 referencetime=me.epoch('UTC','2008/11/24/00:00:00')); 00113 sm.observe(sourcename="fake",spwname='LBand', starttime='-4.0h', stoptime='+4.0h'); 00114 sm.close(); 00115 00116 00117 if(choice in ['XY_all','XY_par','XX_only','YY_only']): 00118 tb.open(basename+'.ms/FEED',nomodify=False); 00119 ff = tb.getcol('POLARIZATION_TYPE'); 00120 ff[0].fill('X'); 00121 ff[1].fill('Y'); 00122 tb.putcol('POLARIZATION_TYPE',ff); 00123 tb.close(); 00124 00125 os.system('rm -rf '+clname); 00126 ###cl.open(); 00127 refRA = qa.unit(ra0); 00128 refDEC = qa.unit(dec0); 00129 cl.addcomponent(flux=stokesvals,fluxunit="Jy", 00130 dir=me.direction(rf='J2000', 00131 v0=qa.add(refRA,"0.0arcmin"),v1=qa.add(refDEC,"0.0arcmin")), 00132 shape="point",freq='1.5GHz', 00133 spectrumtype="spectral index",index=0.0); 00134 cl.rename(filename=clname); 00135 cl.close(); 00136 00137 print >>logfile, 'Predicting : ', clname , ' onto ' , msname; 00138 ft(vis=msname,complist=clname,usescratch=True); 00139 00140 tb.open(msname,nomodify=False); 00141 moddata = tb.getcol(columnname='MODEL_DATA'); 00142 tb.putcol(columnname='DATA',value=moddata); 00143 tb.putcol(columnname='CORRECTED_DATA',value=moddata); 00144 moddata.fill(0.0); 00145 tb.putcol(columnname='MODEL_DATA',value=moddata); 00146 tb.close(); 00147 00148 00149 ##################################### 00150 00151 00152 ### Run clean 00153 def runstokes(msname="Point/point_linXY.ms", imname='sctest',stokes='IQ'): 00154 cmd = 'rm -rf ' + imname + '*'; 00155 os.system(cmd); 00156 print >>logfile, '### Running clean with stokes : ', stokes, ' on MS : ', msname; 00157 clean( vis = msname, 00158 imagename = imname, 00159 niter = 100, 00160 gain = 0.5, 00161 psfmode = "clark", 00162 imagermode = '', 00163 multiscale = [], 00164 imsize = [100], 00165 cell = ['10.0arcsec', '10.0arcsec'], 00166 stokes = stokes, 00167 weighting = "natural" 00168 ); 00169 return checkoutput(msname,imname,stokes); 00170 ################################ 00171 00172 ### Check output images 00173 def checkoutput(msname='Point/point_linXY.ms', imname='sctest', stokes='IQ'): 00174 ## check output stokes coordinates. 00175 resultvals={'ms':msname, 'userstokes':stokes}; 00176 if( not os.path.exists( imname+'.image' ) ): 00177 print >>logfile, "### No output file found"; 00178 resultvals['answer']='no output image'; 00179 return resultvals; 00180 ia.open(imname+'.image'); 00181 csys = ia.coordsys(); 00182 outstokes = csys.stokes(); 00183 if( type(outstokes)==str ): 00184 outstokes = [outstokes]; 00185 midpix = (ia.shape())[0]/2; 00186 rvals={}; 00187 for st in range(0,len(outstokes)): 00188 print >>logfile, '### Output ', outstokes[st], ":", (ia.pixelvalue([midpix,midpix,st,0]))['value']; 00189 rvals[outstokes[st]]=(ia.pixelvalue([midpix,midpix,st,0]))['value']; 00190 ia.close(); 00191 resultvals['answer']=rvals; 00192 return resultvals; 00193 ################################### 00194 00195 ################################### 00196 00197 ########################## 00198 ### Convert IQUV to all stokes dictionary 00199 ########################## 00200 def convertToStokes(stokesvals=[]): 00201 slist = {}; 00202 slist['I'] = stokesvals[0]; 00203 slist['Q']= stokesvals[1]; 00204 slist['U'] = stokesvals[2]; 00205 slist['V']= stokesvals[3]; 00206 slist['RR']= ( slist['I'] + slist['V'] ); 00207 slist['LL']= ( slist['I'] - slist['V'] ); 00208 slist['RL']= ( slist['Q'] + complex(0,slist['U']) ); 00209 slist['LR']= ( slist['Q'] - complex(0,slist['U']) ); 00210 slist['XX']= ( slist['I'] + slist['Q'] ); 00211 slist['YY']= ( slist['I'] - slist['Q'] ); 00212 slist['XY']= ( slist['U'] + complex(0,slist['V']) ); 00213 slist['YX']= ( slist['U'] - complex(0,slist['V']) ); 00214 return slist; 00215 00216 ############################## 00217 00218 00219 ### Run a bunch of tests -- try all combinations allowed at the clean task level. 00220 #### I, IV, IQ, QU, UV, IQUV, RR, LL, RRLL, XX,YY, XXYY 00221 def doAllChecks(stokesvals=[],parentpath='.'): 00222 soptions = ['I','Q','U','V', 'IV', 'QU', 'IQ', 'UV', 'IQU', 'IUV', 'IQUV', 'RR', 'LL', 'RL','LR', 'RRLL','RLLR','XX','YY','XY','YX','XXYY','XYYX']; 00223 choicelist = ['RL_all','RL_par','RR_only','LL_only','XY_all','XY_par','XX_only','YY_only']; 00224 00225 msnames = {'RL_all':'Point/point_linRL.ms' , 'RL_par':'Point/point_parRL.ms' , 'RR_only':'Point/point_onlyRR.ms' ,'LL_only':'Point/point_onlyLL.ms' ,'XY_all':'Point/point_linXY.ms' ,'XY_par':'Point/point_parXY.ms' ,'XX_only':'Point/point_onlyXX.ms' ,'YY_only':'Point/point_onlyYY.ms' } 00226 00227 ## Select what tests to run. 00228 #soptions = ['I','IV','IQU','IQUV']; 00229 #choicelist=['XY_all']; 00230 00231 resultlist=[]; 00232 for mss in choicelist: 00233 for sto in soptions: 00234 print >>logfile, '### ---------------------------------------------------------------- '; 00235 if( not os.path.exists( parentpath+'/'+msnames[mss] ) ): 00236 print >>logfile, '### Did not find MS. Making ', parentpath+'/'+msnames[mss]; 00237 makeMS(mss, stokesvals,parentpath); 00238 else: 00239 print >>logfile, '### Found MS ', parentpath+'/'+msnames[mss]; 00240 resultlist.append( runstokes(msname=parentpath+'/'+msnames[mss],stokes=sto) ); 00241 return resultlist; 00242 ############################################### 00243 00244 ##################################################### 00245 #### END : Helper function definition #### 00246 ##################################################### 00247 00248 ## Now, run the tests. 00249 00250 # Data : Simulated in the test directory. 00251 pathname=os.environ.get('CASAPATH').split()[0] 00252 00253 # Initialize status flag 00254 regstate = True; 00255 00256 # Start timers 00257 startTime=time.time() 00258 startProc=time.clock() 00259 00260 # Run All Tests 00261 parentpath='.'; 00262 stokesvals = [1.0,2.0,3.0,4.0]; 00263 resultlist = doAllChecks(stokesvals=stokesvals,parentpath=parentpath); 00264 00265 # Stop timers 00266 endProc=time.clock() 00267 endTime=time.time() 00268 00269 00270 # Perform the checks 00271 print >>logfile,'' 00272 print >>logfile,'*********************** Comparison of results **********************************' 00273 00274 truestokes=convertToStokes(stokesvals); 00275 00276 countpass=0; 00277 countfail=0; 00278 countwrongnumbers=0; 00279 totalcount=0; 00280 for stok in resultlist: 00281 slist = stok['answer']; 00282 if(slist == 'no output image'): 00283 print >>logfile, 'FAIL : ',stok; 00284 countfail=countfail+1; 00285 else: 00286 stat=True; 00287 for sout in slist: 00288 if( abs( slist[sout]['value'] - truestokes[sout] ) > 1e-03 ): 00289 stat=False; 00290 if(stat==True): 00291 print >>logfile, 'PASS : ', stok; 00292 countpass=countpass+1; 00293 else: 00294 print >>logfile, 'FAIL : WRONG NUMBERS : ', stok; 00295 countfail=countfail+1; 00296 countwrongnumbers=countwrongnumbers+1; 00297 totalcount=totalcount+1; 00298 ################################ 00299 00300 # Truth for initial set of tests 00301 # Feb 18 2011 : This tests clark clean (allows only npol=1,2,4), and stokes=[1.0,2.0,3.0,4.0] 00302 true_totalcount = 184; 00303 true_countpass = 40; 00304 true_countfail = 144; 00305 true_countwrongnumbers = 14; 00306 00307 print >>logfile, '*********************************************************************************'; 00308 print >>logfile, 'OUTPUT -- Total count : ', totalcount , ' --- Passed : ', countpass , ' --- Failed : ', countfail , ' (', countwrongnumbers, ' wrong numbers) '; 00309 print >>logfile, 'TRUTH -- Total count : ', true_totalcount , ' --- Passed : ', true_countpass , ' --- Failed : ', true_countfail , ' (', true_countwrongnumbers, ' wrong numbers) '; 00310 print >>logfile, '*********************************************************************************'; 00311 00312 if( not (totalcount==true_totalcount) or not (countpass==true_countpass) or not (countfail==true_countfail) or not (countwrongnumbers==true_countwrongnumbers) ): 00313 regstate=False; 00314 ################################ 00315 00316 # Final verdict 00317 if(regstate): 00318 print >>logfile,'PASSED regression test for stokes selection in imager' 00319 print '' 00320 print 'Regression PASSED' 00321 print '' 00322 else: 00323 print >>logfile,'FAILED regression test for stokes selection in imager' 00324 print '' 00325 print 'Regression FAILED' 00326 print '' 00327 00328 print >>logfile,'' 00329 00330 # Print timing info 00331 print >>logfile,'********************************************************************************' 00332 print >>logfile,'** Benchmarking **' 00333 print >>logfile,'********************************************************************************' 00334 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00335 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00336 print >>logfile,'Processing rate MB/s was: '+str(278./(endTime - startTime)) 00337 print >>logfile,'* *' 00338 print >>logfile,'********************************************************************************' 00339 00340 logfile.close() 00341 00342
1.8.0
