|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################################################## 00002 # # 00003 # Test Name: # 00004 # imagepoltest_regression.py # 00005 # Rationale for Inclusion: # 00006 # This is a Python translation of the Glish assay test # 00007 # imagepolservertest.g It tests many imagepol tool methods. # 00008 # The first test is a general test of everything at a basic level. # 00009 # Succeeding tests work on individual areas of the imagepol tool. # 00010 # Some of the tests included forced errors. As long as imagepoltest # 00011 # finally returns True it has succeeded (regardless of what error # 00012 # messages you might see). # 00013 # Features Tested: # 00014 # open, done, summary, stokesi, stokesq, stokesu, stokesv, stokes, # 00015 # linpolint, linpolposang, totpolint, fraclinpol, fractotpol, pol, # 00016 # sigmalinpolint, sigmalinpolposang, sigmatotpolint, stokes # 00017 # sigmastokesi, sigmastokesq, sigmastokesu, sigmastokesv, sigmastokes # 00018 # sigmafraclinpol, sigmafractotpol, imagepoltestimage, rotationmeasure # 00019 # fourierrotationmeasure, complexlinpol, complexfraclinpol, makecomplex, # 00020 # depolratio, sigmadepolratio # 00021 # Success/failure criteria: # 00022 # Internally tests each method for pass/fail. Throws an uncaught # 00023 # exception ending test when an unexpected failure occurs. # 00024 # All tests pass if script runs to completion. # 00025 # # 00026 ############################################################################## 00027 # # 00028 # Converted by RRusk 2007-11-16 from imagepolservertest.g # 00029 # # 00030 ############################################################################## 00031 00032 import os 00033 import shutil 00034 import time 00035 import math 00036 import random 00037 from regression_utility import note 00038 00039 po = casac.imagepol() 00040 00041 def imagepoltest(which=None): 00042 # 00043 dowait = true 00044 # 00045 def info(message): 00046 note(message, origin="imagepoltest") 00047 def stop(message): 00048 note(message, priority="SEVERE", origin="imagepoltest") 00049 raise RuntimeError, message 00050 def fail(message=""): 00051 stop(message) 00052 def cleanup(dir): 00053 if (os.path.isdir(dir)): 00054 info("Cleaning up directory "+dir) 00055 def errFunc(raiser, problemPath, excInfo): 00056 note(raiser.__name__+'failed on'+problemPath,"SEVERE") 00057 raise RuntimeError, "Cleanup of " + dir + " fails!" 00058 shutil.rmtree(dir,0,errFunc) 00059 return True 00060 00061 def addnoise(data, sigma): 00062 if len(data.shape)==3: 00063 for i in range(data.shape[0]): 00064 for j in range(data.shape[1]): 00065 for k in range(data.shape[2]): 00066 data[i][j][k] += random.normalvariate(0.0, sigma) 00067 else: 00068 stop('unhandled array shape '+str(data.shape)+' in addnoise') 00069 return True 00070 00071 # 3D only 00072 def make_data(imshape, stokes): 00073 if imshape[2]>4: fail() 00074 if imshape[2]!=len(stokes): fail() 00075 data = ia.makearray(0, [imshape[0], imshape[1], imshape[2]]) 00076 for k in range(imshape[2]): 00077 for i in range(imshape[0]): 00078 for j in range(imshape[1]): 00079 data[i][j][k] = stokes[k] 00080 return data 00081 00082 def alleqnum(x,num,tolerance=0): 00083 if len(x.shape)==3: 00084 for i in range(x.shape[0]): 00085 for j in range(x.shape[1]): 00086 for k in range(x.shape[2]): 00087 if not (abs(x[i][j][k]-num) < tolerance): 00088 print "x[",i,"][",j,"][",k,"]=", x[i][j][k] 00089 return false 00090 else: 00091 stop('unhandled array shape in alleqnum') 00092 return true 00093 00094 def mean(x): 00095 result=0 00096 if len(x.shape)==2: 00097 for i in range(x.shape[0]): 00098 for j in range(x.shape[1]): 00099 result += x[i][j] 00100 result = result/float(x.shape[0]*x.shape[1]) 00101 else: 00102 if len(x.shape)==3: 00103 for i in range(x.shape[0]): 00104 for j in range(x.shape[1]): 00105 for k in range(x.shape[2]): 00106 result += x[i][j][k] 00107 result = result /float(x.shape[0]*x.shape[1]*x.shape[2]) 00108 else: 00109 stop('unhandled array shape '+str(x.shape)+' in mean') 00110 return result 00111 00112 def max_with_location(x): 00113 max=-1e16 00114 loc=[-1,-1,-1,-1] 00115 if len(x.shape)==4: 00116 for i in range(x.shape[0]): 00117 for j in range(x.shape[1]): 00118 for k in range(x.shape[2]): 00119 for l in range(x.shape[3]): 00120 if x[i][j][k][l] > max: 00121 max = x[i][j][k][l] 00122 loc = [i,j,k,l] 00123 else: 00124 stop('unhandled array shape '+str(x.shape)+' in max_with_location') 00125 return max, loc 00126 00127 def test1(): 00128 info('') 00129 info('') 00130 info('') 00131 info('Test 1 - open, done, and summary') 00132 testdir = 'imagepoltest_temp' 00133 if not cleanup(testdir): return False 00134 # Make the directory 00135 try: 00136 os.mkdir(testdir) 00137 except IOError, e: 00138 note(e, "SEVERE") 00139 raise RuntimeError, "mkdir " + testdir + " fails!" 00140 # 00141 # Make RA/DEC only 00142 # 00143 imname = testdir+'/imagefromshape.image' 00144 myim = ia.newimagefromshape(imname, [10,10]) 00145 if not myim: fail() 00146 try: 00147 note("Expect SEVERE error and failure here") 00148 po.open(imname) 00149 isfail = False 00150 except Exception, e: 00151 note("Caught expected Exception") 00152 isfail = True 00153 if not isfail: 00154 stop('open 1 unexpectedly did not fail') 00155 if not myim.remove(True): fail() 00156 # 00157 # RA/DEC/I 00158 # 00159 myim = ia.newimagefromshape(imname, [10,10,1]) 00160 if not myim: fail() 00161 try: 00162 note("Expect SEVERE error and failure here") 00163 po.open(imname) 00164 isfail = False 00165 except Exception, e: 00166 note("Caught expected Exception") 00167 isfail = True 00168 if not isfail: 00169 stop('open 2 unexpectedly did not fail') 00170 if not myim.remove(True): fail() 00171 # 00172 # RA/DEC/IQ 00173 # 00174 myim = ia.newimagefromshape(imname, [10,10,2]) 00175 if not myim: fail() 00176 try: 00177 note("Expect SEVERE error and failure here") 00178 po.open(imname) 00179 isfail = False 00180 except Exception, e: 00181 note("Caught expected Exception") 00182 isfail = True 00183 if not isfail: 00184 stop('open 3 unexpectedly did not fail') 00185 if not myim.remove(True): fail() 00186 # 00187 # RA/DEC/IQU 00188 # 00189 myim = ia.newimagefromshape(imname, [10,10,3]) 00190 if not myim: fail() 00191 try: 00192 po.open(imname) 00193 isfail = False 00194 except Exception, e: 00195 isfail = True 00196 if isfail: stop('open 4 failed') 00197 if not po.done(): fail() 00198 if not myim.remove(True): fail() 00199 # 00200 # RA/DEC/IQUV 00201 # 00202 myim = ia.newimagefromshape(imname, [10,10,4]) 00203 if not myim: fail() 00204 try: 00205 po.open(imname) 00206 isfail = False 00207 except Exception, e: 00208 isfail = True 00209 if isfail: stop('open 5 failed') 00210 if not po.done(): fail() 00211 # 00212 # Test tool constructor 00213 # 00214 if not(po.open(myim.torecord())): 00215 return stop('open 6 failed') 00216 if not po.done(): fail() 00217 if not(po.open(myim.name())): 00218 return stop('open 7 failed') 00219 00220 if not myim.done(): fail() 00221 # 00222 # Utility functions 00223 # 00224 ok = po.summary() 00225 if not ok: fail() 00226 if not po.done(): fail() 00227 # 00228 return cleanup(testdir) 00229 00230 00231 def test2(): 00232 info('') 00233 info('') 00234 info('') 00235 info('Test 2 - stokesi, stokesq, stokesu, stokesv, stokes') 00236 testdir = 'imagepoltest_temp' 00237 if not cleanup(testdir): return False 00238 # Make the directory 00239 try: 00240 os.mkdir(testdir) 00241 except IOError, e: 00242 note(e, "SEVERE") 00243 raise RuntimeError, "mkdir " + testdir + " fails!" 00244 00245 # 00246 # Make some data 00247 # 00248 shape = [10,10,4] 00249 stokes = [14,2,3,4] 00250 data = make_data(shape, stokes) 00251 if len(data)==0: fail() 00252 00253 # 00254 # Make image - RA/DEC/IQUV 00255 # 00256 imname = testdir+'/imagefromarray.image' 00257 myim = ia.newimagefromarray(imname, data) 00258 if not myim: fail() 00259 if not myim.done(): fail() 00260 00261 # 00262 # Attach image to polarimetry tool 00263 # 00264 if not po.open(imname): fail() 00265 00266 # 00267 # Get Stokes images 00268 # 00269 s = po.stokesi() 00270 if not s: fail() 00271 pixels = s.getchunk() 00272 if not alleqnum(pixels, stokes[0], 0.0001): 00273 stop('Stokes I values are wrong') 00274 if not s.done(): fail() 00275 # 00276 s = po.stokesq() 00277 if not s: fail() 00278 pixels = s.getchunk() 00279 if not alleqnum(pixels, stokes[1], 0.0001): 00280 stop('Stokes Q values are wrong') 00281 if not s.done(): fail() 00282 # 00283 s = po.stokesu() 00284 if not s: fail() 00285 pixels = s.getchunk() 00286 if not alleqnum(pixels, stokes[2], 0.0001): 00287 stop('Stokes U values are wrong') 00288 if not s.done(): fail() 00289 # 00290 s = po.stokesv(); 00291 if not s: fail() 00292 pixels = s.getchunk() 00293 if not alleqnum(pixels, stokes[3], 0.0001): 00294 stop('Stokes V values are wrong') 00295 if not s.done(): fail() 00296 00297 ss = ["i","q","u","v"] 00298 for i in range(len(ss)): 00299 s = po.stokes(ss[i]) 00300 if not s: fail() 00301 pixels = s.getchunk() 00302 if not alleqnum(pixels, stokes[i], 0.0001): 00303 stop('Stokes '+ss[i]+' values are wrong') 00304 if not s.done(): fail() 00305 00306 # 00307 try: 00308 note('Expect SEVERE error and Exception here') 00309 s = po.stokes('fish') 00310 isfail = False 00311 except Exception, e: 00312 note('Cuaght expected Exception') 00313 isfail = True 00314 if not isfail: 00315 stop('Function stokes unexpectedly did not fail') 00316 00317 # 00318 if not po.done(): fail() 00319 # 00320 return cleanup(testdir) 00321 00322 00323 def test3(): 00324 info('') 00325 info('') 00326 info('') 00327 info('Test 3 - linpolint, linpolposang, totpolint') 00328 # 00329 testdir = 'imagepoltest_temp' 00330 if not cleanup(testdir): return False 00331 # Make the directory 00332 try: 00333 os.mkdir(testdir) 00334 except IOError, e: 00335 note(e, "SEVERE") 00336 raise RuntimeError, "mkdir " + testdir + " fails!" 00337 00338 # 00339 # Make some data 00340 # 00341 shape = [256,256,4] 00342 stokes = [14,2,3,4] 00343 data = make_data(shape, stokes) 00344 sigma = 0.01 * stokes[1] 00345 ok = addnoise(data, sigma) 00346 if not ok: fail() 00347 00348 # 00349 # Make image - RA/DEC/IQUV 00350 # 00351 imname = testdir+'/imagefromarray.image' 00352 myim = ia.newimagefromarray(imname, data) 00353 if not myim: fail() 00354 if not myim.done(): fail() 00355 00356 # 00357 # Open image with polarimetry tool 00358 # 00359 if not po.open(imname): fail() 00360 00361 # 00362 # Linearly polarized intensity 00363 # 00364 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + sigma*sigma) 00365 s = po.linpolint(debias=F) 00366 if not s: fail() 00367 pixels = s.getchunk() 00368 d = abs(mean(pixels)-pp) / float(pp) 00369 if d > 0.01: 00370 stop('Linearly polarized intensity values are wrong') 00371 if not s.done(): fail() 00372 # 00373 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] - sigma*sigma) 00374 s = po.linpolint(debias=T, clip=10.0) 00375 if not s: fail() 00376 pixels = s.getchunk() 00377 d = abs(mean(pixels)-pp) / float(pp) 00378 if d > 0.01: 00379 stop('Debiased linearly polarized intensity values (1) are wrong') 00380 if not s.done(): fail() 00381 # 00382 s = po.linpolint(debias=T, clip=10.0, sigma=sigma) 00383 if not s: fail() 00384 pixels = s.getchunk() 00385 d = abs(mean(pixels)-pp) / float(pp) 00386 if d > 0.01: 00387 stop('Debiased linearly polarized intensity values (2) are wrong') 00388 if not s.done(): fail() 00389 00390 # 00391 # Linearly polarized position angle 00392 pp = (180.0/(2.0*math.pi)) * math.atan2(stokes[2], stokes[1]) 00393 s = po.linpolposang() 00394 if not s: fail() 00395 pixels = s.getchunk() 00396 d = abs(mean(pixels)-pp) / float(pp) 00397 if d > 0.01: 00398 stop('Linearly polarized position angles are wrong') 00399 if not s.done(): fail() 00400 00401 # 00402 # 00403 # Total polarized intensity 00404 # 00405 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + 00406 stokes[3]*stokes[3] + sigma*sigma) 00407 s = po.totpolint(debias=F) 00408 if not s: fail() 00409 pixels = s.getchunk() 00410 d = abs(mean(pixels)-pp) / float(pp) 00411 if d>0.1: 00412 stop('Total polarized intensity values are wrong') 00413 if not s.done(): fail() 00414 # 00415 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + 00416 stokes[3]*stokes[3] - sigma*sigma) 00417 s = po.totpolint(debias=T, clip=10.0) 00418 if not s: fail() 00419 pixels = s.getchunk() 00420 d = abs(mean(pixels)-pp) / float(pp) 00421 if d>0.1: 00422 stop('Debiased total polarized intensity values (1) are wrong') 00423 if not s.done(): fail() 00424 # 00425 s = po.totpolint(debias=T, clip=10.0, sigma=sigma) 00426 if not s: fail() 00427 pixels = s.getchunk() 00428 d = abs(mean(pixels)-pp) / float(pp) 00429 if d>0.1: 00430 stop('Debiased total polarized intensity values (2) are wrong') 00431 if not s.done(): fail() 00432 # 00433 if not po.done(): fail() 00434 # 00435 return cleanup(testdir) 00436 00437 def test4(): 00438 info('') 00439 info('') 00440 info('') 00441 info('Test 4 - fraclinpol, fractotpol') 00442 # 00443 testdir = 'imagepoltest_temp' 00444 if not cleanup(testdir): return False 00445 # Make the directory 00446 try: 00447 os.mkdir(testdir) 00448 except IOError, e: 00449 note(e, "SEVERE") 00450 raise RuntimeError, "mkdir " + testdir + " fails!" 00451 # 00452 # Make some data 00453 # 00454 shape = [256,256,4] 00455 stokes = [14,2,3,4] 00456 data = make_data(shape, stokes) 00457 sigma = 0.01 * stokes[1] 00458 ok = addnoise(data, sigma) 00459 if not ok: fail() 00460 00461 # 00462 # Make image - RA/DEC/IQUV 00463 # 00464 imname = testdir+'/imagefromarray.image' 00465 myim = ia.newimagefromarray(imname, data) 00466 if not myim: fail() 00467 if not myim.done(): fail() 00468 00469 # 00470 # Open image with polarimetry tool 00471 # 00472 if not po.open(imname): fail() 00473 00474 # 00475 # Fractional linear polarization 00476 # 00477 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + 00478 sigma*sigma) / float(stokes[0]) 00479 s = po.fraclinpol(debias=F) 00480 if not s: fail() 00481 pixels = s.getchunk() 00482 d = abs(mean(pixels)-pp) / float(pp) 00483 if d>0.1: 00484 stop('Fractional linear polarization values are wrong') 00485 if not s.done(): fail() 00486 # 00487 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] - 00488 sigma*sigma) / float(stokes[0]) 00489 s = po.fraclinpol(debias=T, clip=10.0) 00490 if not s: fail() 00491 pixels = s.getchunk() 00492 d = abs(mean(pixels)-pp) / float(pp) 00493 if d>0.1: 00494 stop('Debiased fractional linear polarization values (1) are wrong') 00495 if not s.done(): fail() 00496 # 00497 s = po.fraclinpol(debias=T, clip=10.0, sigma=sigma) 00498 if not s: fail() 00499 pixels = s.getchunk() 00500 d = abs(mean(pixels)-pp) / float(pp) 00501 if d>0.1: 00502 stop('Debiased fractional linear polarization values (2) are wrong') 00503 if not s.done(): fail() 00504 00505 # 00506 # Fractional total polarization 00507 # 00508 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + 00509 stokes[3]*stokes[3] + sigma*sigma) / float(stokes[0]) 00510 s = po.fractotpol(debias=F) 00511 if not s: fail() 00512 pixels = s.getchunk() 00513 d = abs(mean(pixels)-pp) / float(pp) 00514 if d>0.1: 00515 stop('Fractional total polarization values are wrong') 00516 if not s.done(): fail() 00517 # 00518 pp = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + 00519 stokes[3]*stokes[3] - sigma*sigma) / float(stokes[0]) 00520 s = po.fractotpol(debias=T, clip=10.0) 00521 if not s: fail() 00522 pixels = s.getchunk() 00523 d = abs(mean(pixels)-pp) / float(pp) 00524 if d>0.1: 00525 stop('Debiased fractional total polarization values (1) are wrong') 00526 if not s.done(): fail() 00527 # 00528 s = po.fractotpol(debias=T, clip=10.0, sigma=sigma) 00529 if not s: fail() 00530 pixels = s.getchunk() 00531 d = abs(mean(pixels)-pp) / float(pp) 00532 if d>0.1: 00533 stop('Debiased fractional total polarization values (2) are wrong') 00534 if not s.done(): fail() 00535 # 00536 if not po.done(): fail() 00537 # 00538 return cleanup(testdir) 00539 00540 def test5(): 00541 info('') 00542 info('') 00543 info('') 00544 info('Test 5 - pol') 00545 # 00546 testdir = 'imagepoltest_temp' 00547 if not cleanup(testdir): return False 00548 # Make the directory 00549 try: 00550 os.mkdir(testdir) 00551 except IOError, e: 00552 note(e, "SEVERE") 00553 raise RuntimeError, "mkdir " + testdir + " fails!" 00554 # 00555 # Make some data 00556 # 00557 shape = [256,256,4] 00558 stokes = [14,2,3,4] 00559 data = make_data(shape, stokes) 00560 sigma = 0.01 * stokes[1] 00561 ok = addnoise(data, sigma) 00562 if not ok: fail() 00563 00564 # 00565 # Make image - RA/DEC/IQUV 00566 # 00567 imname = testdir+'/imagefromarray.image' 00568 myim = ia.newimagefromarray(imname, data) 00569 if not myim: fail() 00570 if not myim.done(): fail() 00571 00572 # 00573 # Open image with polarimetry tool 00574 # 00575 if not po.open(imname): fail() 00576 00577 # 00578 # We just test that each function runs, not its results, as this 00579 # just packages previously tested functions 00580 # 00581 which = ["lpi","tpi","lppa","flp","ftp"] 00582 for i in which: 00583 s = po.pol(i, debias=F) 00584 if not s: fail() 00585 if not s.done(): fail() 00586 # 00587 s = po.pol(i, debias=T, clip=10.0) 00588 if not s: fail() 00589 if not s.done(): fail() 00590 # 00591 s = po.pol(i, debias=T, clip=10.0, sigma=sigma) 00592 if not s: fail() 00593 if not s.done(): fail() 00594 try: 00595 note("Expect SEVERE error and Exception here") 00596 s = po.pol('fish') 00597 isfail = False 00598 except Exception, e: 00599 note("Caught expected Exception") 00600 isfail = True 00601 if not isfail: 00602 stop('Function pol unexpectedly did not fail') 00603 00604 # 00605 if not po.done(): fail() 00606 return cleanup(testdir) 00607 00608 def test6(): 00609 info('') 00610 info('') 00611 info('') 00612 info('Test 6 - sigmalinpolint, sigmalinpolposang, sigmatotpolint') 00613 # 00614 testdir = 'imagepoltest_temp' 00615 if not cleanup(testdir): return False 00616 # Make the directory 00617 try: 00618 os.mkdir(testdir) 00619 except IOError, e: 00620 note(e, "SEVERE") 00621 raise RuntimeError, "mkdir " + testdir + " fails!" 00622 # 00623 # Make some data 00624 # 00625 shape = [256,256,4] 00626 stokes = [14,2,3,4] 00627 data = make_data(shape, stokes) 00628 sigma = 0.01 * stokes[1] 00629 ok = addnoise(data, sigma) 00630 if not ok: fail() 00631 00632 # 00633 # Make image - RA/DEC/IQUV 00634 # 00635 imname = testdir+'/imagefromarray.image' 00636 myim = ia.newimagefromarray(imname, data) 00637 if not myim: fail() 00638 if not myim.done(): fail() 00639 00640 # 00641 # Open image with polarimetry tool 00642 # 00643 if not po.open(imname): fail() 00644 00645 # 00646 # Error in linearly polarized intensity 00647 # 00648 s = po.sigmalinpolint(clip=10.0) 00649 if not s: fail() 00650 d = abs(s-sigma)/float(sigma) 00651 if d>0.01: 00652 stop('Sigma for linearly polarized intensity (1) is wrong') 00653 # 00654 s = po.sigmalinpolint(clip=10.0, sigma=sigma) 00655 if not s: fail() 00656 d = abs(s-sigma)/float(sigma) 00657 if d>0.1: 00658 stop('Sigma for linearly polarized intensity (2) is wrong') 00659 00660 # 00661 # Error in linearly polarized position angle 00662 # 00663 s = po.sigmalinpolposang(clip=10.0) 00664 if not s: fail() 00665 data = s.getchunk() 00666 if not s.done(): fail() 00667 lpi = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2]) 00668 s2 = 180.0 * sigma / float(lpi) / 2.0 / math.pi 00669 if not abs((mean(data)-s2)/float(s2)) < 0.01: 00670 stop('Sigma for linearly polarized position angle (1) is wrong') 00671 # 00672 s = po.sigmalinpolposang(clip=10.0, sigma=sigma) 00673 if not s: fail() 00674 data = s.getchunk() 00675 if not s.done(): fail() 00676 lpi = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2]) 00677 s2 = 180.0 * sigma / float(lpi) / 2.0 / math.pi 00678 d = abs(mean(data)-s2)/float(s2) 00679 if d>0.1: 00680 stop('Sigma for linearly polarized position angle (2) is wrong') 00681 00682 # 00683 # Error in total linearly polarized intensity 00684 # 00685 s = po.sigmatotpolint(clip=10.0) 00686 if not s: fail() 00687 d = abs(s-sigma)/float(sigma) 00688 if d>0.1: 00689 stop('Sigma for total polarized intensity (1) is wrong') 00690 # 00691 s = po.sigmatotpolint(clip=10.0, sigma=sigma) 00692 if not s: fail() 00693 d = abs(s-sigma)/float(sigma) 00694 if d>0.1: 00695 stop('Sigma for total polarized intensity (2) is wrong') 00696 00697 # 00698 # Cleanup 00699 # 00700 if not po.done(): fail() 00701 return cleanup(testdir) 00702 00703 def test7(): 00704 info('') 00705 info('') 00706 info('') 00707 info('Test 7 - sigma, sigmastokesi, sigmastokesq, sigmastokesu, sigmastokesv, sigmastokes') 00708 # 00709 testdir = 'imagepoltest_temp' 00710 if not cleanup(testdir): return False 00711 # Make the directory 00712 try: 00713 os.mkdir(testdir) 00714 except IOError, e: 00715 note(e, "SEVERE") 00716 raise RuntimeError, "mkdir " + testdir + " fails!" 00717 # 00718 # Make some data 00719 # 00720 shape = [256,256,4] 00721 stokes = [14,2,3,4] 00722 data = make_data(shape, stokes) 00723 sigma = 0.01 * stokes[1] 00724 ok = addnoise(data, sigma) 00725 if not ok: fail() 00726 00727 # 00728 # Make image - RA/DEC/IQUV 00729 # 00730 imname = testdir+'/imagefromarray.image' 00731 myim = ia.newimagefromarray(imname, data) 00732 if not myim: fail() 00733 if not myim.done(): fail() 00734 00735 # 00736 # Open image with polarimetry tool 00737 # 00738 if not po.open(imname): fail() 00739 00740 # 00741 # Best guess at thermal noise 00742 # 00743 s = po.sigma(clip=100.0) 00744 if not s: fail() 00745 d = abs(s-sigma)/float(sigma) 00746 if d>0.1: 00747 stop('Sigma is wrong') 00748 00749 # 00750 # Error in stokes I, Q, U, V 00751 # 00752 s = po.sigmastokesi(clip=100.0) 00753 if not s: fail() 00754 d = abs(s-sigma)/float(sigma) 00755 if d>0.1: 00756 stop('Sigma for Stokes I is wrong') 00757 # 00758 s = po.sigmastokesq(clip=100.0) 00759 if not s: fail() 00760 if not abs(s-sigma) < 0.001: 00761 stop('Sigma for Stokes Q is wrong') 00762 # 00763 s = po.sigmastokesu(clip=100.0) 00764 if not s: fail() 00765 d = abs(s-sigma)/float(sigma) 00766 if d>0.1: 00767 stop('Sigma for Stokes U is wrong') 00768 # 00769 s = po.sigmastokesv(clip=100.0) 00770 if not s: fail() 00771 d = abs(s-sigma)/float(sigma) 00772 if d>0.1: 00773 stop('Sigma for Stokes V is wrong') 00774 # 00775 which = ["I","Q","U","V"] 00776 for i in which: 00777 s = po.sigmastokes(which=i, clip=100.0) 00778 if not s: fail() 00779 try: 00780 note("Expect SEVERE error and Exception here") 00781 s = po.sigmastokes(which='fish') 00782 isfail = false 00783 except Exception, e: 00784 note("Caught expected exception") 00785 isfail = true 00786 if not isfail: 00787 stop('Function sigmastokes unexpectedly did not fail') 00788 00789 # 00790 # Cleanup 00791 # 00792 if not po.done(): fail() 00793 return cleanup(testdir) 00794 00795 def test8(): 00796 info('') 00797 info('') 00798 info('') 00799 info('Test 8 - sigmafraclinpol, sigmafractotpol') 00800 # 00801 testdir = 'imagepoltest_temp' 00802 if not cleanup(testdir): return False 00803 # Make the directory 00804 try: 00805 os.mkdir(testdir) 00806 except IOError, e: 00807 note(e, "SEVERE") 00808 raise RuntimeError, "mkdir " + testdir + " fails!" 00809 # 00810 # Make some data 00811 # 00812 shape = [256,256,4] 00813 stokes = [14,2,3,4] 00814 data = make_data(shape, stokes) 00815 sigma = 0.01 * stokes[1] 00816 ok = addnoise(data, sigma) 00817 if not ok: fail() 00818 00819 # 00820 # Make image - RA/DEC/IQUV 00821 # 00822 imname = testdir+'/imagefromarray.image' 00823 myim = ia.newimagefromarray(imname, data) 00824 if not myim: fail() 00825 if not myim.done(): fail() 00826 00827 # 00828 # Open image with polarimetry tool 00829 # 00830 if not po.open(imname): fail() 00831 00832 # 00833 # Error in fractional linearly polarized intensity 00834 # 00835 s = po.sigmafraclinpol(clip=10.0) 00836 if not s: fail() 00837 data = s.getchunk() 00838 if not s.done(): fail() 00839 pi = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2]) 00840 m = pi / float(stokes[0]) 00841 s2 = m * math.sqrt( (sigma/float(pi))*(sigma/float(pi)) + 00842 (sigma/float(stokes[0]))*(sigma/float(stokes[0]))) 00843 d = abs(mean(data)-s2)/float(s2) 00844 if d>0.1: 00845 stop('Sigma for fractional linear polarization (1) is wrong') 00846 # 00847 s = po.sigmafraclinpol(clip=10.0, sigma=sigma) 00848 if not s: fail() 00849 data = s.getchunk() 00850 if not s.done(): fail() 00851 pi = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2]) 00852 m = pi / float(stokes[0]) 00853 s2 = m * math.sqrt( (sigma/float(pi))*(sigma/float(pi)) + 00854 (sigma/float(stokes[0]))*(sigma/float(stokes[0]))) 00855 d = abs(mean(data)-s2)/float(s2) 00856 if d>0.1: 00857 stop('Sigma for fractional linear polarization (2) is wrong') 00858 00859 # 00860 # Error in fractional total polarized intensity 00861 # 00862 s = po.sigmafractotpol(clip=10.0) 00863 if not s: fail() 00864 data = s.getchunk() 00865 if not s.done(): fail() 00866 pi = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + 00867 stokes[3]*stokes[3]) 00868 m = pi / float(stokes[0]) 00869 s2 = m * math.sqrt( (sigma/float(pi))*(sigma/float(pi)) + 00870 (sigma/float(stokes[0]))*(sigma/float(stokes[0]))) 00871 d = abs(mean(data)-s2)/float(s2) 00872 if d>0.1: 00873 stop('Sigma for fractional total polarization (1) is wrong') 00874 # 00875 s = po.sigmafractotpol(clip=10.0, sigma=sigma) 00876 if not s: fail() 00877 data = s.getchunk() 00878 if not s.done(): fail() 00879 pi = math.sqrt(stokes[1]*stokes[1] + stokes[2]*stokes[2] + 00880 stokes[3]*stokes[3]) 00881 m = pi / float(stokes[0]) 00882 s2 = m * math.sqrt( (sigma/float(pi))*(sigma/float(pi)) + 00883 (sigma/float(stokes[0]))*(sigma/float(stokes[0]))) 00884 d = abs(mean(data)-s2)/float(s2) 00885 if d>0.1: 00886 stop('Sigma for fractional total polarization (2) is wrong') 00887 00888 # 00889 # Cleanup 00890 # 00891 if not po.done(): fail() 00892 return cleanup(testdir) 00893 00894 def test9(): 00895 info('') 00896 info('') 00897 info('') 00898 info('Test 9 - imagepoltestimage, rotationmeasure') 00899 # 00900 testdir = 'imagepoltest_temp' 00901 if not cleanup(testdir): return False 00902 # Make the directory 00903 try: 00904 os.mkdir(testdir) 00905 except IOError, e: 00906 note(e, "SEVERE") 00907 raise RuntimeError, "mkdir " + testdir + " fails!" 00908 00909 # 00910 # Make imagepoltestimage 00911 # 00912 outfile = testdir+'/imagepoltestimage.image' 00913 rm = 200.0 00914 pa0 = 10.0 00915 sigma = 0.001 00916 nf = 32 00917 f0 = 1.4e9 00918 bw = 128.0e6 00919 if not po.imagepoltestimage (outfile=outfile, rm=rm, pa0=pa0, 00920 sigma=sigma, nx=32, ny=32, nf=nf, 00921 f0=f0, bw=bw): fail() 00922 # 00923 # Rotation Measure 00924 # 00925 rmmax = rm + 100.0 00926 maxpaerr = 100000.0 00927 rmfg = 0.0 00928 rmname = testdir+'/rm.image' 00929 rmename = testdir+'/rme.image' 00930 pa0name = testdir+'/pa0.image' 00931 pa0ename = testdir+'/pa0e.image' 00932 ok = po.rotationmeasure(rm=rmname, pa0=pa0name, 00933 rmerr=rmename, pa0err=pa0ename, 00934 sigma=sigma, rmfg=rmfg, 00935 rmmax=rmmax, maxpaerr=maxpaerr) 00936 if not ok: fail() 00937 00938 # 00939 # Check results 00940 # 00941 rmim = ia.newimagefromfile(rmname) 00942 rmeim = ia.newimagefromfile(rmename) 00943 pa0im = ia.newimagefromfile(pa0name) 00944 pa0eim = ia.newimagefromfile(pa0ename) 00945 # 00946 err = mean(rmeim.getchunk()) 00947 diff = mean(rmim.getchunk()) - rm 00948 if abs(diff) > 3*err: 00949 stop('Recovered RM is wrong') 00950 # 00951 err = mean(pa0eim.getchunk()) 00952 diff = mean(pa0im.getchunk()) - pa0 00953 if abs(diff) > 3*err: 00954 stop('Recovered Position Angle at zero frequency is wrong') 00955 00956 # 00957 # Cleanup 00958 # 00959 if not po.done(): fail() 00960 if not rmim.done(): fail() 00961 if not rmeim.done(): fail() 00962 if not pa0im.done(): fail() 00963 if not pa0eim.done(): fail() 00964 # 00965 return cleanup(testdir) 00966 00967 def test10(): 00968 info('') 00969 info('') 00970 info('') 00971 info('Test 10 - fourierrotationmeasure') 00972 # 00973 testdir = 'imagepoltest_temp' 00974 if not cleanup(testdir): return False 00975 # Make the directory 00976 try: 00977 os.mkdir(testdir) 00978 except IOError, e: 00979 note(e, "SEVERE") 00980 raise RuntimeError, "mkdir " + testdir + " fails!" 00981 00982 # 00983 # Make imagepoltestimage 00984 # 00985 outfile = testdir+'/imagepoltestimage.image' 00986 rm = 1.0e6 00987 pa0 = 0.0 00988 sigma = 0.0001 00989 nf = 512 00990 f0 = 1.4e9 00991 bw = 15.625e3 * nf 00992 if not po.imagepoltestimage (outfile=outfile, rm=rm, pa0=pa0, 00993 sigma=sigma, nx=1, ny=1, nf=nf, 00994 f0=f0, bw=bw): fail() 00995 # 00996 # Fourier Rotation Measure 00997 # 00998 ampname = testdir+'/amp.image' 00999 ok = po.fourierrotationmeasure(amp=ampname) 01000 if not ok: fail() 01001 01002 # 01003 # Check results 01004 # 01005 ampim = ia.newimagefromfile(ampname) 01006 srec = ampim.summary(list=F) 01007 rminc = srec['incr'][3] 01008 rmrefpix = srec['refpix'][3] 01009 idx = int((rm+rminc/2.0)/float(rminc) + rmrefpix) 01010 y = ampim.getchunk() 01011 max,loc = max_with_location(y) 01012 if idx != loc[3]: 01013 stop('Peak of RM spectrum is in wrong channel') 01014 01015 # 01016 # Cleanup 01017 # 01018 if not po.done(): fail() 01019 if not ampim.done(): fail() 01020 # 01021 return cleanup(testdir) 01022 01023 def test11(): 01024 info('') 01025 info('') 01026 info('') 01027 info('Test 11 - complexlinpol, complexfraclinpol, makecomplex') 01028 # 01029 testdir = 'imagepoltest_temp' 01030 if not cleanup(testdir): return False 01031 # Make the directory 01032 try: 01033 os.mkdir(testdir) 01034 except IOError, e: 01035 note(e, "SEVERE") 01036 raise RuntimeError, "mkdir " + testdir + " fails!" 01037 # 01038 # Make some data 01039 # 01040 shape = [2,2,4] 01041 stokes = [14,2,3,4] 01042 data = make_data(shape, stokes) 01043 01044 # 01045 # Make image - RA/DEC/IQUV 01046 # 01047 imname = testdir+'/imagefromarray.image' 01048 myim = ia.newimagefromarray(imname, data) 01049 if not myim: fail() 01050 if not myim.done(): fail() 01051 01052 # 01053 # Open image with polarimetry tool 01054 # 01055 if not po.open(imname): fail() 01056 01057 # 01058 # Complex linear polarization 01059 # 01060 s = testdir+'/complexpol.image' 01061 ok = po.complexlinpol(s) 01062 if not ok: fail() 01063 # 01064 expr = 'real("'+s+'")' 01065 rim = ia.imagecalc(pixels=expr) 01066 if not rim: fail() 01067 qpixels = rim.getchunk() 01068 if len(qpixels)==0: fail() 01069 if not rim.done(): fail() 01070 d = abs(mean(qpixels)-stokes[1]) 01071 if d > 0.001: 01072 stop('Complex linear polarization (1) values (Q) are wrong') 01073 # 01074 expr = 'imag("'+s+'")' 01075 iim = ia.imagecalc(pixels=expr) 01076 if not iim: fail() 01077 upixels = iim.getchunk() 01078 if len(upixels)==0: fail() 01079 if not iim.done(): fail() 01080 d = abs(mean(upixels)-stokes[2]) 01081 if d > 0.001: 01082 stop('Complex linear polarization (1) values (U) are wrong') 01083 01084 # 01085 # Complex fractional polarization 01086 # 01087 s = testdir+'/complexfracpol.image' 01088 ok = po.complexfraclinpol(s) 01089 if not ok: fail() 01090 # 01091 expr = 'real("'+s+'")' 01092 rim = ia.imagecalc(pixels=expr) 01093 if not rim: fail() 01094 qpixels = rim.getchunk() 01095 if len(qpixels)==0: fail() 01096 if not rim.done(): fail() 01097 d = abs(mean(qpixels)-(float(stokes[1])/float(stokes[0]))) 01098 if d > 0.001: 01099 stop('Complex fractional polarization (1) values (Q) are wrong') 01100 # 01101 expr = 'imag("'+s+'")' 01102 iim = ia.imagecalc(pixels=expr) 01103 if not iim: fail() 01104 upixels = iim.getchunk() 01105 if len(upixels) ==0: fail() 01106 if not iim.done(): fail() 01107 d = abs(mean(upixels)-(stokes[2]/float(stokes[0]))) 01108 if d > 0.001: 01109 stop('Complex fractional polarization (1) values (U) are wrong') 01110 01111 # 01112 # Makecomplex 01113 # 01114 q = po.stokesq() 01115 qs = testdir+'/q.image' 01116 q2 = q.subimage(qs) 01117 if not q.done(): fail() 01118 if not q2.done(): fail() 01119 # 01120 u = po.stokesu() 01121 us = testdir+'/u.image' 01122 u2 = u.subimage(us) 01123 if not u.done(): fail() 01124 if not u2.done(): fail() 01125 # 01126 lpi = po.linpolint() 01127 lpis = testdir+'/lpi.image' 01128 lpi2 = lpi.subimage(lpis) 01129 if not lpi.done(): fail() 01130 if not lpi2.done(): fail() 01131 # 01132 lppa = po.linpolposang() 01133 lppas = testdir+'/lppa.image' 01134 lppa2 = lppa.subimage(lppas) 01135 if not lppa.done(): fail() 01136 if not lppa2.done(): fail() 01137 # 01138 s = testdir+'/cplx1.image' 01139 po.makecomplex(s, real=qs, imag=us) 01140 # 01141 expr = 'real("'+s+'")' 01142 rim = ia.imagecalc(pixels=expr) 01143 if not rim: fail() 01144 rpixels = rim.getchunk() 01145 if len(rpixels)==0: fail() 01146 if not rim.done(): fail() 01147 d = abs(mean(rpixels)-(stokes[1])) 01148 if d > 0.001: 01149 stop('Complex linear polarization (2) values (Q) are wrong') 01150 # 01151 expr = 'imag("'+s+'")' 01152 iim = ia.imagecalc(pixels=expr) 01153 if not iim: fail() 01154 ipixels = iim.getchunk() 01155 if len(ipixels)==0: fail() 01156 if not iim.done(): fail() 01157 d = abs(mean(ipixels)-(stokes[2])) 01158 if d > 0.001: 01159 stop('Complex linear polarization (2) values (U) are wrong') 01160 # 01161 s = testdir+'/cplx2.image' 01162 po.makecomplex(s, amp=lpis, phase=lppas) 01163 # 01164 expr = 'real("'+s+'")' 01165 rim = ia.imagecalc(pixels=expr) 01166 if not rim: fail() 01167 rpixels = rim.getchunk() 01168 if len(rpixels)==0: fail() 01169 if not rim.done(): fail() 01170 d = abs(mean(rpixels)-(stokes[1])) 01171 if d > 0.001: 01172 stop('Complex linear polarization (3) values (Q) are wrong') 01173 # 01174 expr = 'imag("'+s+'")' 01175 iim = ia.imagecalc(pixels=expr) 01176 if not iim: fail() 01177 ipixels = iim.getchunk() 01178 if len(ipixels)==0: fail() 01179 if not iim.done(): fail() 01180 d = abs(mean(ipixels)-(stokes[2])) 01181 if d > 0.001: 01182 stop('Complex linear polarization (3) values (U) are wrong') 01183 # 01184 if not po.done(): fail() 01185 return cleanup(testdir) 01186 01187 def test12(): 01188 info('') 01189 info('') 01190 info('') 01191 info('Test 12 - depolratio, sigmadepolratio') 01192 # 01193 testdir = 'imagepoltest_temp' 01194 if not cleanup(testdir): return False 01195 # Make the directory 01196 try: 01197 os.mkdir(testdir) 01198 except IOError, e: 01199 note(e, "SEVERE") 01200 raise RuntimeError, "mkdir " + testdir + " fails!" 01201 # 01202 # Make some data 01203 # 01204 shape = [256,256,4] 01205 # 01206 stokes1 = [14,2,3,0.2] 01207 data1 = make_data(shape, stokes1) 01208 sigma1 = 0.01 * stokes1[1] 01209 if not addnoise(data1, sigma1): fail() 01210 # 01211 stokes2 = [13,1.5,2.8,0.1] 01212 data2 = make_data(shape, stokes2) 01213 sigma2 = 0.01 * stokes2[1] 01214 if not addnoise(data2, sigma2): fail() 01215 01216 # 01217 # Make images - RA/DEC/IQUV 01218 # 01219 imname1 = testdir+'/imagefromarray.image1' 01220 myim1 = ia.newimagefromarray(imname1, data1) 01221 if not myim1: fail() 01222 if not myim1.done(): fail() 01223 # 01224 imname2 = testdir+'/imagefromarray.image2' 01225 myim2 = ia.newimagefromarray(imname2, data2) 01226 if not myim2: fail() 01227 if not myim2.done(): fail() 01228 # 01229 # Make polarimetry server 01230 # 01231 if not po.open(imname1): fail() 01232 01233 # 01234 # Depolarization ratio 01235 # 01236 i1 = stokes1[0] 01237 i2 = stokes1[0] 01238 ei1 = sigma1 01239 ei2 = sigma2 01240 # 01241 p1 = math.sqrt(stokes1[1]*stokes1[1] + stokes1[2]*stokes1[2]) 01242 p2 = math.sqrt(stokes2[1]*stokes2[1] + stokes2[2]*stokes2[2]) 01243 ep1 = sigma1 01244 ep2 = sigma2 01245 # 01246 m1 = p1 / float(stokes1[0]) 01247 m2 = p2 / float(stokes2[0]) 01248 em1 = m1 * math.sqrt( (ep1*ep1/float(p1*p1)) + (ei1*ei1/float(i1*i1))) 01249 em2 = m2 * math.sqrt( (ep2*ep2/float(p2*p2)) + (ei2*ei2/float(i2*i2))) 01250 # 01251 dd = m1 / float(m2) 01252 edd = dd * math.sqrt( (em1*em1/float(m1*m1)) + (em2*em2/float(m2*m2)) ) 01253 # 01254 depol = po.depolratio(infile=imname2, debias=F); # Use file name 01255 if not depol: fail() 01256 pixels = depol.getchunk() 01257 diff = abs(mean(pixels)-dd) / float(dd) 01258 if diff > 0.01: 01259 stop('Depolarization ratio values are wrong') 01260 if not depol.done(): fail() 01261 # 01262 myim2 = ia.newimagefromfile(imname2) 01263 if not myim2: fail() 01264 depol = po.depolratio(infile=myim2.name(), debias=F); # Use Image tool 01265 if not depol: fail() 01266 pixels = depol.getchunk() 01267 diff = abs(mean(pixels)-dd) / float(dd) 01268 if diff > 0.01: 01269 stop('Depolarization ratio values are wrong') 01270 if not myim2.done(): fail() 01271 if not depol.done(): fail() 01272 # 01273 # Error in depolarization ratio 01274 # 01275 edepol = po.sigmadepolratio(infile=imname2, debias=F); # Use file name 01276 if not edepol: fail() 01277 pixels = edepol.getchunk() 01278 diff = abs(mean(pixels)-edd) / float(edd) 01279 if diff > 0.01: 01280 stop('Depolarization ratio error values are wrong') 01281 if not edepol.done(): fail() 01282 01283 # 01284 if not po.done(): fail() 01285 return cleanup(testdir) 01286 01287 #run tests 01288 test1() 01289 test2() 01290 test3() 01291 test4() 01292 test5() 01293 test6() 01294 test7() 01295 test8() 01296 test9() 01297 test10() 01298 test11() 01299 test12() 01300 print '' 01301 print 'Regression PASSED' 01302 print '' 01303 01304 01305 01306 ##################################################################### 01307 # 01308 # Get on with it 01309 # 01310 note ('', priority='WARN', origin='imagepoltest_regression.py') 01311 note ('These tests include intentional errors. You should expect to see error messages in the log.', 01312 priority='WARN', origin='imagepoltest.py') 01313 note ('If the test finally returns True, then it has succeeded\n\n', 01314 priority='WARN', origin='imagepoltest.py') 01315 note ('', priority='WARN', origin='imagepoltest.py') 01316 # 01317 01318 Benchmarking = True 01319 if Benchmarking: 01320 startTime = time.time() 01321 regstate = False 01322 for i in range(1): 01323 imagepoltest() 01324 endTime = time.time() 01325 regstate = True 01326 else: 01327 imagepoltest() 01328 01329 #exit()
1.8.0
