casa: regressions/m51_3sim_regression.py Source File

Go to the documentation of this file.
00001 ###############################################
00002 # Regression Script for simdata of a 2d image #
00003 # (single dish only simulation)               #
00004 ###############################################
00005 
00006 import os, time
00007 
00008 #modelname="M51HA.MODEL"
00009 modelname="m51ha.model"
00010 if os.path.exists(modelname):
00011     shutil.rmtree(modelname)
00012 
00013 noise=False # add noise
00014 
00015 projname = "m51c"
00016 
00017 startTime = time.time()
00018 startProc = time.clock()
00019 
00020 print '--Running simlation of M51 (ALMA-12m INT + ACA-7m INT + 12m TP) --'
00021 # configs are in the repository
00022 
00023 l=locals() 
00024 if not l.has_key("repodir"): 
00025     repodir=os.getenv("CASAPATH").split(' ')[0]
00026 
00027 print casa['build']
00028 print 'I think the data repository is at '+repodir
00029 datadir=repodir+"/data/regression/simdata/"
00030 cfgdir=repodir+"/data/alma/simmos/"
00031 #importfits(fitsimage=datadir+modelname,imagename="m51.image")
00032 shutil.copytree(datadir+modelname,modelname)
00033 
00034 
00035 
00036 #======================================
00037 # 12m INT
00038 
00039 default("simobserve")
00040 
00041 project = projname
00042 # Clear out results from previous runs.
00043 #os.system('rm -rf '+project+'*')
00044 os.system('rm -rf '+project)
00045 
00046 skymodel = modelname
00047 inbright = '0.004'
00048 indirection = 'B1950 23h59m59.96 -34d59m59.50'
00049 incell = '0.1arcsec'
00050 incenter = '330.076GHz'
00051 inwidth = '50MHz'
00052 
00053 setpointings = True
00054 integration = '10s'
00055 mapsize = '1arcmin'
00056 maptype = "hex"
00057 pointingspacing = '9arcsec'
00058 
00059 #observe = True
00060 obsmode = "int"
00061 refdate='2012/11/21/20:00:00'
00062 totaltime = '3600s'
00063 #sdantlist = cfgdir+'aca.tp.cfg'
00064 #sdant = 0
00065 
00066 antennalist="alma;0.5arcsec"
00067 
00068 if noise:
00069     thermalnoise = 'tsys-atm'  #w/ noise 
00070     user_pwv=3.0
00071 
00072 if not l.has_key('interactive'): interactive=False
00073 if interactive:
00074     graphics="both"
00075 else:
00076     graphics="file"
00077 
00078 verbose=True
00079 overwrite=True
00080 
00081 inp()
00082 go()
00083 
00084 
00085 
00086 
00087 #========================================================
00088 # 12m TP
00089 
00090 
00091 default("simobserve")
00092 
00093 project = projname
00094 
00095 skymodel = modelname
00096 inbright = '0.004'
00097 indirection = 'B1950 23h59m59.96 -34d59m59.50'
00098 incell = '0.1arcsec'
00099 incenter = '330.076GHz'
00100 inwidth = '50MHz'
00101 
00102 setpointings = True
00103 integration = '10s'
00104 mapsize = '1arcmin'
00105 #maptype = "hex"
00106 maptype = "square"
00107 pointingspacing = '9arcsec'
00108 
00109 #observe = True
00110 obsmode = "sd"
00111 refdate='2012/11/21/20:00:00'
00112 totaltime = '2h'
00113 sdantlist = cfgdir+'aca.tp.cfg'
00114 sdant = 0
00115 #antennalist=""
00116 
00117 if noise:
00118     thermalnoise = 'tsys-atm'  #w/ noise 
00119     user_pwv=3.0
00120 
00121 if not l.has_key('interactive'): interactive=False
00122 if interactive:
00123     graphics="both"
00124 else:
00125     graphics="file"
00126 
00127 verbose=True
00128 overwrite=True
00129 
00130 inp()
00131 go()
00132 
00133 
00134 
00135 
00136 #========================================================
00137 # ACA
00138 
00139 
00140 default("simobserve")
00141 
00142 project = projname
00143 
00144 skymodel = modelname
00145 inbright = '0.004'
00146 indirection = 'B1950 23h59m59.96 -34d59m59.50'
00147 incell = '0.1arcsec'
00148 incenter = '330.076GHz'
00149 inwidth = '50MHz'
00150 
00151 setpointings = True
00152 integration = '10s'
00153 mapsize = '1arcmin'
00154 maptype = "hex"
00155 #maptype = "square"
00156 pointingspacing = '15arcsec'
00157 
00158 #observe = True
00159 obsmode = "int"
00160 refdate='2012/11/21/20:00:00'
00161 totaltime = '3' # times through the map
00162 #sdantlist = cfgdir+'aca.tp.cfg'
00163 #sdant = 0
00164 
00165 antennalist="aca.i.cfg"
00166 
00167 if noise:
00168     thermalnoise = 'tsys-atm'  #w/ noise 
00169     user_pwv=3.0
00170 
00171 if not l.has_key('interactive'): interactive=False
00172 if interactive:
00173     graphics="both"
00174 else:
00175     graphics="file"
00176 
00177 verbose=True
00178 overwrite=True
00179 
00180 inp()
00181 go()
00182 
00183 
00184 
00185 
00186 #==============================================
00187 # clean
00188 
00189 
00190 default("simanalyze")
00191 
00192 project = projname
00193 
00194 # clean ACA with SD model
00195 
00196 image = True
00197 if noise:
00198     vis = '$project.aca.i.noisy.ms,$project.aca.tp.sd.noisy.ms'  #w/ noise
00199 else:
00200     vis = '$project.aca.i.ms,$project.aca.tp.sd.ms'  #w/ noise
00201 imsize = [512,512]
00202 #imdirection = 'B1950 23h59m59.96 -34d59m59.50'
00203 cell = '0.2arcsec'
00204 
00205 analyze = True
00206 # show psf & residual are not available for SD-only simulation
00207 showpsf = False
00208 showresidual = False
00209 showconvolved = True
00210 
00211 if not l.has_key('interactive'): interactive=False
00212 if interactive:
00213     graphics="both"
00214 else:
00215     graphics="file"
00216 
00217 inp()
00218 go()
00219 
00220 
00221 
00222 
00223 default("simanalyze")
00224 
00225 project = projname
00226 
00227 # clean ALMA with ACA+SD model
00228 
00229 image = True
00230 if noise:
00231     vis = '$project.alma_0.5arcsec.noisy.ms'
00232 else:
00233     vis = '$project.alma_0.5arcsec.ms'
00234 imsize = [512,512]
00235 #imdirection = 'B1950 23h59m59.96 -34d59m59.50'
00236 cell = '0.2arcsec'
00237 modelimage="$project.aca.i.image"
00238 
00239 analyze = True
00240 # show psf & residual are not available for SD-only simulation
00241 showpsf = False
00242 showresidual = False
00243 showconvolved = True
00244 
00245 if not l.has_key('interactive'): interactive=False
00246 if interactive:
00247     graphics="both"
00248 else:
00249     graphics="file"
00250 
00251 inp()
00252 go()
00253 
00254 
00255 
00256 
00257 
00258 
00259 endTime = time.time()
00260 endProc = time.clock()
00261 
00262 # Regression
00263 
00264 test_name = """simdata observation of M51 (ALMA-12m INT + ACA-7m INT + 12m TP)"""
00265 
00266 ia.open(project+"/"+project + '.alma_0.5arcsec.image')
00267 m51both_stats=ia.statistics(verbose=False,list=False)
00268 ia.close()
00269 
00270 
00271 ia.open(project+"/"+project + '.alma_0.5arcsec.diff')
00272 m51both_diffstats=ia.statistics(verbose=False,list=False)
00273 ia.close()
00274 
00275 # # reference statistic values for simulated image
00276 # refstats = { 'max': 0.12078,
00277 #              'min': -0.022069,
00278 #              'rms': 0.016695,
00279 #              'sigma': 0.014399,
00280 #              'sum': 951.08 }
00281 
00282 # reference statistic values for simulated image
00283 refstats = { 'max': 0.13314,
00284              'min': -0.023702,
00285              'rms': 0.020402,
00286              'sigma': 0.016707,
00287              'sum': 1402.5 }
00288 
00289 
00290 # # reference statistic values for diff image
00291 # diffstats = {'max': 0.030648,
00292 #              'min': -0.075413,
00293 #              'rms': 0.0096487,
00294 #              'sigma': 0.0096155,
00295 #              'sum': -90.078 }
00296 
00297 # reference statistic values for diff image
00298 diffstats = {'max': 0.03206,
00299              'min': -0.10785,
00300              'rms': 0.014134,
00301              'sigma': 0.013433,
00302              'sum': -526.38 }
00303 
00304 # relative tolerances to reference values
00305 reftol   = {'sum':  1e-2,
00306             'max':  1e-2,
00307             'min':  1e-2,
00308             'rms':  1e-2,
00309             'sigma': 1e-2}
00310 
00311 import datetime
00312 datestring = datetime.datetime.isoformat(datetime.datetime.today())
00313 outfile    = project+"/"+project + '.' + datestring + '.log'
00314 logfile    = open(outfile, 'w')
00315 
00316 print 'Writing regression output to ' + outfile + "\n"
00317 print >> logfile,casa['build']
00318 
00319 loghdr = """
00320 ********** Regression *****************
00321 """
00322 
00323 print >> logfile, loghdr
00324 
00325 # more info
00326 ms.open(project+"/"+project+".alma_0.5arcsec.ms")
00327 print >> logfile, "Noiseless MS, amp stats:"
00328 print >> logfile, ms.statistics('DATA','amp')
00329 print >> logfile, "Noiseless MS, phase stats:"
00330 print >> logfile, ms.statistics('DATA','phase')
00331 ms.close()
00332 #ms.open(project+".noisy.ms")
00333 #print >> logfile, "Noisy MS, amp stats:"
00334 #print >> logfile, ms.statistics('DATA','amp')
00335 #print >> logfile, "Noisy MS, phase stats:"
00336 #print >> logfile, ms.statistics('DATA','phase')
00337 #ms.close()
00338 
00339 
00340 regstate = True
00341 rskes = refstats.keys()
00342 rskes.sort()
00343 for ke in rskes:
00344     adiff=abs(m51both_stats[ke][0] - refstats[ke])/abs(refstats[ke])
00345     if adiff < reftol[ke]:
00346         print >> logfile, "* Passed %-5s image test, got % -11.5g expected % -11.5g." % (ke, m51both_stats[ke][0], refstats[ke])
00347     else:
00348         print >> logfile, "* FAILED %-5s image test, got % -11.5g instead of % -11.5g." % (ke, m51both_stats[ke][0], refstats[ke])
00349         regstate = False
00350 
00351 rskes = diffstats.keys()
00352 rskes.sort()
00353 for ke in rskes:
00354     adiff=abs(m51both_diffstats[ke][0] - diffstats[ke])/abs(diffstats[ke])
00355     if adiff < reftol[ke]:
00356         print >> logfile, "* Passed %-5s  diff test, got % -11.5g expected % -11.5g." % (ke, m51both_diffstats[ke][0], diffstats[ke])
00357     else:
00358         print >> logfile, "* FAILED %-5s  diff test, got % -11.5g instead of % -11.5g." % (ke, m51both_diffstats[ke][0], diffstats[ke])
00359         regstate = False
00360 
00361 # this script doesn't have sensible values yet 20100928
00362 regstate=True        
00363 
00364 print >> logfile,'---'
00365 if regstate:
00366     print >> logfile, 'Passed',
00367     print ''
00368     print 'Regression PASSED'
00369     print ''
00370 else:
00371     print >> logfile, 'FAILED',
00372     print ''
00373     print 'Regression FAILED'
00374     print ''
00375 
00376 print >> logfile, 'regression test for simdata of M51 (ALMA-12m INT + ACA-7m INT + 12m TP).'
00377 print >>logfile,'---'
00378 print >>logfile,'*********************************'
00379     
00380 print >>logfile,''
00381 print >>logfile,'********** Benchmarking **************'
00382 print >>logfile,''
00383 print >>logfile,'Total wall clock time was: %8.3f s.' % (endTime - startTime)
00384 print >>logfile,'Total CPU        time was: %8.3f s.' % (endProc - startProc)
00385 print >>logfile,'Wall processing  rate was: %8.3f MB/s.' % (17896.0 /
00386                                                             (endTime - startTime))
00387 
00388 ### Get last modification time of .ms.
00389 msfstat = os.stat(project+"/"+project+'.alma_0.5arcsec.ms')
00390 print >>logfile,'* Breakdown:                           *'
00391 print >>logfile,'*  generating visibilities took %8.3fs,' % (msfstat[8] - startTime)
00392 print >>logfile,'*************************************'
00393     
00394 logfile.close()
00395                                                     
00396 print '--Finished simdata of M51 (total power+interferometric) regression--'