|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################### 00002 # # 00003 # Regression/Benchmarking Script for NGC 4826 # 00004 # Single Dish + Synthesis # 00005 ############################################### 00006 00007 import time 00008 import os 00009 00010 os.system('rm -rf n4826_t* mosaic.flux mosaic.model mosaic.image mosaic.residual n12m_gaussian.im') 00011 00012 pathname=os.environ.get('CASAPATH').split()[0] 00013 datapath = pathname + '/data/regression/ATST3/NGC4826/' 00014 00015 print '--Copy data to local directory--' 00016 mspath='cp -r '+datapath+'n4826_both.ms .' 00017 os.system(mspath) 00018 os.system('chmod -R a+wx n4826_both.ms') 00019 00020 startTime = time.time() 00021 startProc = time.clock() 00022 00023 00024 print '--Feather--' 00025 # Starting from: 00026 # BIMA mosaic image (Moment 0) : n4826_mom0.im 00027 # NRAO 12m OTF image (Moment 0): n4826_12mmom0.im 00028 default('feather') 00029 feather('n4826_tfeather.im',datapath+'n4826_mom0.im',datapath+'n4826_12mmom0.im') 00030 feathertime = time.time() 00031 #combo: Max:1.533498e+02 Flux:1.523515e+03 Jy rms:1.187669e+01 00032 #Pcombo:Max:1.511604e+02 Flux:2.016790e+03 Jy rms:1.873132e+01 00033 #BIMA: Max:1.627327e+02 Flux:1.035314e+04 Jy rms:1.587101e+01 00034 #12m: Max:9.939910e+02 Flux:2.600573e+03 Jy rms:1.159703e+02 00035 00036 print '--Feather cube--' 00037 # Starting from: 00038 # BIMA mosaic image (Moment 0) : n4826_bima.im 00039 # NRAO 12m OTF image (Moment 0): NGC4826.12motf.chan.fits 00040 default('importfits') 00041 importfits(datapath+'NGC4826.12motf.chan.fits','n4826_t12mchan.im') 00042 default('imhead') 00043 imhead('n4826_t12mchan.im',mode='put',hdkey='bunit',hdvalue='Jy/beam') 00044 #imhead('n4826_t12mchan.im',mode='put',hdkey='beam',hdvalue='55arcsec, 55arcsec, 0deg') 00045 ###wow what an idea...more to type to change beam 00046 imhead('n4826_t12mchan.im',mode='put',hdkey='beammajor',hdvalue='55arcsec') 00047 imhead('n4826_t12mchan.im',mode='put',hdkey='beamminor',hdvalue='55arcsec') 00048 imhead('n4826_t12mchan.im',mode='put',hdkey='beampa',hdvalue='0deg') 00049 00050 00051 default('feather') 00052 feather('n4826_tfeathercube.im',datapath+'n4826_bima.im','n4826_t12mchan.im') 00053 feathercubetime = time.time() 00054 # Maximum value 1.881626e+00 at [113, 115, 1, 25] (12:56:45.388, +21.40.51.100, I, 1.15172e+11Hz) 00055 # Flux density = 1.054628e+02 Jy 00056 00057 print '--Feather cube - create synth image--' 00058 # Starting from: 00059 # BIMA calibrated visibilities: n4826_both.ms 00060 # NRAO 12m OTF cube: NGC4826.12motf.chan.fits 00061 default('clean') 00062 clean(vis='n4826_both.ms', imagename='mosaic', 00063 nchan=30, start=45, width=4, spw='0~2', 00064 field='0~6', 00065 cell=[1., 1.], imsize=[256, 256], 00066 stokes='I', 00067 mode='channel', 00068 interpolation='nearest', 00069 psfmode='clark', 00070 imagermode='mosaic', 00071 ftmachine='mosaic', 00072 niter=500, 00073 gain=0.1, # Specified for certainty, 6/10/2009 00074 scaletype='SAULT', minpb=0.1) # Changed from 0.01 6/10/2009 00075 default('importfits') 00076 importfits(datapath+'NGC4826.12motf.chan.fits', 'n4826_t12mchan2.im') 00077 default('imhead') 00078 imhead('n4826_t12mchan2.im', mode='put', hdkey='bunit', hdvalue='Jy/beam') 00079 #imhead('n4826_t12mchan2.im',mode='put',hdkey='beam',hdvalue='55arcsec,55arcsec,0deg') 00080 imhead('n4826_t12mchan2.im',mode='put',hdkey='beammajor',hdvalue='55arcsec') 00081 imhead('n4826_t12mchan2.im',mode='put',hdkey='beamminor',hdvalue='55arcsec') 00082 imhead('n4826_t12mchan2.im',mode='put',hdkey='beampa',hdvalue='0deg') 00083 default('feather') 00084 feather('n4826_tfeathersynth.im','mosaic.image','n4826_t12mchan2.im') 00085 #Combo: Max: 1.881626e+00 Flux: 1.054628e+02 Jy Rms: 8.911050e-02 00086 #Pcombo: Max: 1.681789e+00 Flux: 1.454376e+02 Jy Rms: 8.196454e-02 00087 feathersynthtime = time.time() 00088 00089 print '--Single Dish as Model--' 00090 # Starting from: 00091 # BIMA calibrated visibilities: n4826_both.ms 00092 # NRAO 12m OTF cube: n4826_t12mchan.im 00093 default('clean') 00094 clean(vis='n4826_both.ms', imagename='n4826_tjoint1', 00095 nchan=30, start=45, width=4, spw='0~2', 00096 field='0~6', 00097 cell=[1., 1.], imsize=[256, 256], 00098 stokes='I', 00099 mode='channel', 00100 psfmode='clark', 00101 niter=500, 00102 cyclefactor=4, 00103 interpolation='nearest', 00104 imagermode='mosaic', 00105 ftmachine='mosaic', 00106 modelimage='n4826_t12mchan.im', 00107 scaletype='SAULT', minpb=0.1) 00108 sdmodeltime = time.time() 00109 #Combo: Max: 1.891995e+00 Flux: 1.391578e+02 Jy Rms: 9.056366e-02 00110 #Pcombo: Max: 1.681789e+00 Flux: 1.454376e+02 Jy Rms: 8.196454e-02 00111 # im:=image('n4826_joint1.restored'); 00112 # im_mom0:=im.moments(outfile='n4826_tmom0.im', moments=0, axis=4, 00113 # mask='indexin(4,[5:26])',includepix=[0.070,1000.0]); 00114 #Combo: Max: 1.891995e+00 Flux:1.098632e+04 Jy Rms: 1.678347e+01 00115 #BIMA: Max:1.627327e+02 Flux:1.035314e+04 Jy rms:1.587101e+01 00116 #Pcombo:Max:1.511604e+02 Flux:2.016790e+03 Jy rms:1.873132e+01 00117 00118 print '--Joint Deconvolution--' 00119 # Regrid 12m image onto synth imaging coordinates 00120 # Tool-kit only (currently) 00121 ia.open('n4826_tjoint1.image') 00122 csys=ia.coordsys() 00123 ia.close() 00124 #### regrid SD image so that shape and coordinate parameters match 00125 ia.open('n4826_t12mchan.im') 00126 ia.regrid(outfile='n4826_t12motf.chregrid.im',shape=[256,256,1,30],csys=csys.torecord()) 00127 ia.close() 00128 #### deconvolve SD image with a guess of PB with msclean 00129 dc.open('n4826_t12motf.chregrid.im', psf='') 00130 dc.makegaussian('n12m_gaussian.im' ,bmaj='55arcsec', bmin='55arcsec', bpa='0deg', 00131 normalize=false) 00132 dc.close() 00133 dc.open('n4826_t12motf.chregrid.im', psf='n12m_gaussian.im') 00134 dc.setscales(scalemethod='uservector', uservector=[30., 60.]) 00135 dc.clean(algorithm='msclean', model='n4826_tjoint2', niter=100, gain=0.3) 00136 #dc.close() 00137 default('clean') 00138 ##### Mosaic the interferometer data...use model from obtain from deconvolve 00139 ##### SD image as starting model 00140 clean(vis='n4826_both.ms', imagename='n4826_tjoint2', 00141 nchan=30, start=45, width=4, spw='0~2', 00142 field='0~6', 00143 cell=[1., 1.], imsize=[256, 256], 00144 stokes='I', 00145 mode='channel', 00146 psfmode='clark', 00147 imagermode='mosaic', 00148 ftmachine='mosaic', 00149 interpolation='nearest', 00150 niter=500, 00151 cyclefactor=4, 00152 modelimage='n4826_tjoint2', 00153 scaletype='SAULT', minpb=0.1) 00154 00155 jointtime = time.time() 00156 00157 endProc = time.clock() 00158 endTime = time.time() 00159 00160 # Regression 00161 00162 test_name_feather1 = 'NGC4826--Feather (Synthesis and SD Moment 0 images' 00163 test_name_feather2 = 'NGC4826--Feather Cube (Synthesis and SD data cube' 00164 test_name_feather3 = 'NGC4826--Create Synth cube; feather with SD FITS cube' 00165 test_name_jc1 = 'NGC4826--Joint deconvolution with SD FITS cube as model' 00166 test_name_jc2 = 'NGC4826--Joint deconvolution with deconvolved SD FITS cube' 00167 00168 ia.open('n4826_tfeather.im') 00169 ia.setbrightnessunit('Jy/beam') 00170 ia.close() 00171 00172 feather1_im=ia.open('n4826_tfeather.im') 00173 f1_stats=ia.statistics();ia.close() 00174 feather2_im=ia.open('n4826_tfeathercube.im') 00175 f2_stats=ia.statistics();ia.close() 00176 feather3_im=ia.open('n4826_tfeathersynth.im') 00177 f3_stats=ia.statistics();ia.close() 00178 jcfeather1_im=ia.open('n4826_tjoint1.image') 00179 jc1_stats=ia.statistics();ia.close() 00180 jcfeather2_im=ia.open('n4826_tjoint2.image') 00181 jc2_stats=ia.statistics() 00182 00183 feather1_immax=f1_stats['max'][0] 00184 feather1_flux=f1_stats['flux'][0] 00185 feather2_immax=f2_stats['max'][0] 00186 feather2_flux=f2_stats['flux'][0] 00187 feather3_immax=f3_stats['max'][0] 00188 feather3_flux=f3_stats['flux'][0] 00189 jc1_immax=jc1_stats['max'][0] 00190 joint1_flux=jc1_stats['flux'][0] 00191 jc2_immax=jc2_stats['max'][0] 00192 joint2_flux=jc2_stats['flux'][0] 00193 00194 00195 #Note flux values differ from AIPS++; now correcting for primary beam so raises 00196 #the noise at the edges; should do stats on central region (soon) 00197 ###these numbers are really fragile to minor changes...like the total flux 00198 ### makes no sense at all. better crieria needed 00199 f1_max=153.3498 # < 12/1/2009 00200 f1_flux=1523.515 # < 12/1/2009 00201 f2_max=1.8816 # < 12/1/2009 00202 f2_flux=105.4628 # < 12/1/2009 00203 00204 #f3_max=1.67 00205 #f3_max=1.52 00206 f3_max=1.60 # 12/1/2009. Are we converging? 00207 #f3_max=1.6825 # 3/17/2010 00208 #f3_max = 1.457553 # 3/19/2010, after increasing clean's start by width/2 00209 00210 f3_flux=110.6 # change in channel centres 05/02/2011 00211 #f3_flux=110.44684 # 3/17/2010 00212 #f3_flux=99.4 # 3/19/2010 (adjusted clean's start by +1) 00213 #f3_flux=83.27245 # 3/19/2010 (adjusted clean's start by another +1 (width/2 total)) 00214 00215 jc1_max=1.67 # < 12/1/2009 00216 #jc1_max=1.71 00217 00218 #jc1_flux=168.87 00219 jc1_flux=230.3 # changes due to channel centre change 00220 #jc1_flux=223.033 # 3/17/2010 00221 #jc1_flux=224.495 # 3/19/2010 (adjusted clean's start by +1) 00222 #jc1_flux=212.652 # 3/19/2010 (adjusted clean's start by another +1 (width/2 total) 00223 # Note that it's almost back to 12/1/2009.) 00224 00225 #jc2_max=1.68 00226 jc2_max=1.53 # < 12/1/2009 00227 00228 #jc2_flux=67.27 00229 jc2_flux=147.7 # changes due to change in channel centre 00230 #jc2_flux=144.9498 # 3/17/2010 00231 #jc2_flux=145.09 # 3/19/2010 (adjusted clean's start by +1) 00232 #jc2_flux=135.89555 # 3/19/2010 (adjusted clean's start by another +1 (width/2 total)) 00233 00234 diff_f1=abs((f1_max-feather1_immax)/f1_max) 00235 diff_f1f=abs((f1_flux-feather1_flux)/f1_flux) 00236 diff_f2=abs((f2_max-feather2_immax)/f2_max) 00237 diff_f2f=abs((f2_flux-feather2_flux)/f2_flux) 00238 diff_f3=abs((f3_max-feather3_immax)/f3_max) 00239 diff_f3f=abs((f3_flux-feather3_flux)/f3_flux) 00240 diff_jc1=abs((jc1_max-jc1_immax)/jc1_max) 00241 diff_jc1f=abs((jc1_flux-joint1_flux)/jc1_flux) 00242 diff_jc2=abs((jc2_max-jc2_immax)/jc2_max) 00243 diff_jc2f=abs((jc2_flux-joint2_flux)/jc2_flux) 00244 00245 import datetime 00246 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00247 outfile='ngc4826c.'+datestring+'.log' 00248 logfile=open(outfile, 'w') 00249 00250 print >>logfile, '' 00251 print >>logfile, '********** Data Summary *********' 00252 print >>logfile, '' 00253 print >>logfile, '* Observer: Project: t108c115.n48' 00254 print >>logfile, '*Observation: BIMA(10 antennas)' 00255 print >>logfile, '* Telescope Observation Date Observer Project' 00256 print >>logfile, '* BIMA [ 4.39941e+09, 4.39942e+09] t108c115.n48' 00257 print >>logfile, '*Data records: 109260 Total integration time = 628492 seconds' 00258 print >>logfile, '* Observed from 04:04:31 to 10:39:23' 00259 print >>logfile, '*Fields: 7' 00260 print >>logfile, '* ID Name Right Ascension Declination Epoch' 00261 print >>logfile, '* 0 NGC4826-F0 12:56:44.24 +21.41.05.10 J2000' 00262 print >>logfile, '* 1 NGC4826-F1 12:56:41.08 +21.41.05.10 J2000' 00263 print >>logfile, '* 2 NGC4826-F2 12:56:42.66 +21.41.43.20 J2000' 00264 print >>logfile, '* 3 NGC4826-F3 12:56:45.82 +21.41.43.20 J2000' 00265 print >>logfile, '* 4 NGC4826-F4 12:56:47.40 +21.41.05.10 J2000' 00266 print >>logfile, '* 5 NGC4826-F5 12:56:45.82 +21.40.27.00 J2000' 00267 print >>logfile, '* 6 NGC4826-F6 12:56:42.66 +21.40.27.00 J2000' 00268 print >>logfile, '*Spectral Windows: (4 unique spectral windows and 1 unique polarization setups)' 00269 print >>logfile, '* SpwID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs' 00270 print >>logfile, '* 0 64 LSRD 114950.387 1562.5 100000 115271.2 YY' 00271 print >>logfile, '* 1 64 LSRD 115040.402 1562.5 100000 115271.2 YY' 00272 print >>logfile, '* 2 64 LSRD 115130.143 1562.5 100000 115271.2 YY' 00273 print >>logfile, '* 3 64 LSRD 115220.157 1562.5 100000 115271.2 YY' 00274 print >>logfile, '*Antennas: 10' 00275 print >>logfile, '* ID= 1-3: ANT1=UNKNOWN, ANT2=UNKNOWN, ANT3=UNKNOWN, ANT4=UNKNOWN,' 00276 print >>logfile,'* ID= 5-7: ANT5=UNKNOWN, ANT6=UNKNOWN, ANT7=UNKNOWN, ANT8=UNKNOWN,' 00277 print >>logfile,'* ID= 9-9: ANT9=UNKNOWN, ANT10=UNKNOWN' 00278 print >>logfile,'*********************************' 00279 print >>logfile,'' 00280 print >>logfile,'********** Regression ***********' 00281 print >>logfile,'* *' 00282 status = {False: "! FAILED", True: "* Passed"} 00283 print >>logfile, status[diff_f1 < 0.05], 'Feather 1 image max test ' 00284 print >>logfile,'*-- Feather 1: Image max '+str(feather1_immax)+','+str(f1_max) 00285 print >>logfile, status[diff_f2 < 0.05], 'Feather 2 image max test' 00286 print >>logfile,'*-- Feather 2: Image max '+str(feather2_immax)+','+str(f2_max) 00287 print >>logfile, status[diff_f3 < 0.05], 'Feather 3 image max test' 00288 print >>logfile,'*-- Feather 3: Image max '+str(feather3_immax)+','+str(f3_max) 00289 print >>logfile, status[diff_jc1 < 0.05], 'Joint Deconvolution 1 image max test' 00290 print >>logfile,'*-- Joint Decon1: Image max '+str(jc1_immax)+','+str(jc1_max) 00291 print >>logfile, status[diff_jc2 < 0.05], 'Joint Deconvolution 2 image max test' 00292 print >>logfile,'*-- Joint Decon2: Image max '+str(jc2_immax)+','+str(jc2_max) 00293 print >>logfile, status[diff_f1f < 0.05], 'Feather 1 flux test ' 00294 print >>logfile,'*-- Feather 1: Flux '+str(feather1_flux)+','+str(f1_flux) 00295 print >>logfile, status[diff_f2f < 0.05], 'Feather 2 flux test' 00296 print >>logfile,'*-- Feather 2: Flux '+str(feather2_flux)+','+str(f2_flux) 00297 print >>logfile, status[diff_f3f < 0.05], 'Feather 3 flux test' 00298 print >>logfile,'*-- Feather 3: Flux '+str(feather3_flux)+','+str(f3_flux) 00299 print >>logfile, status[diff_jc1f < 0.05], 'Joint Deconvolution flux test' 00300 print >>logfile,'*-- Joint Decon1: Flux '+str(joint1_flux)+','+str(jc1_flux) 00301 print >>logfile, status[diff_jc2f < 0.05], 'Joint Deconvolution flux test' 00302 print >>logfile,'*-- Joint Decon2: Flux '+str(joint2_flux)+','+str(jc2_flux) 00303 00304 if (diff_f1 < 0.05 and 00305 diff_f2 < 0.05 and 00306 diff_f3 <0.05 and 00307 diff_jc1 < 0.05 and 00308 diff_f1f < 0.05 and 00309 diff_f2f < 0.05 and 00310 diff_f3f < 0.05 and 00311 diff_jc1f < 0.05): 00312 regstate=True 00313 print >>logfile,'---' 00314 print >>logfile,'Passed Regression test for NGC4826' 00315 print >>logfile,'---' 00316 else: 00317 regstate=False 00318 print >>logfile,'----FAILED Regression test for NGC4826' 00319 print >>logfile,'*********************************' 00320 00321 print >>logfile,'' 00322 print >>logfile,'' 00323 print >>logfile,'********* Benchmarking *****************' 00324 print >>logfile,'* *' 00325 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00326 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00327 print >>logfile,'Processing rate MB/s was: '+str(760./(endTime - startTime)) 00328 # n4826_both.ms x 3 plus image sizesish 00329 print >>logfile,'* Breakdown: *' 00330 print >>logfile,'* feather time was: '+str(feathertime-startTime) 00331 print >>logfile,'* feathercube time was: '+str(feathercubetime-feathertime) 00332 print >>logfile,'* feathersynth time was: '+str(feathersynthtime-feathercubetime) 00333 print >>logfile,'* sdmodel time was: '+str(sdmodeltime-feathersynthtime) 00334 print >>logfile,'* joint decon time was: '+str(jointtime-sdmodeltime) 00335 print >>logfile,'*****************************************' 00336 00337 logfile.close()
1.8.0
