|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ########################################################################## 00002 # imfit_test.py 00003 # 00004 # Copyright (C) 2008, 2009 00005 # Associated Universities, Inc. Washington DC, USA. 00006 # 00007 # This script is free software; you can redistribute it and/or modify it 00008 # under the terms of the GNU Library General Public License as published by 00009 # the Free Software Foundation; either version 2 of the License, or (at your 00010 # option) any later version. 00011 # 00012 # This library is distributed in the hope that it will be useful, but WITHOUT 00013 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00014 # FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00015 # License for more details. 00016 # 00017 # You should have received a copy of the GNU Library General Public License 00018 # along with this library; if not, write to the Free Software Foundation, 00019 # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00020 # 00021 # Correspondence concerning AIPS++ should be adressed as follows: 00022 # Internet email: aips2-request@nrao.edu. 00023 # Postal address: AIPS++ Project Office 00024 # National Radio Astronomy Observatory 00025 # 520 Edgemont Road 00026 # Charlottesville, VA 22903-2475 USA 00027 # 00028 # <author> 00029 # Dave Mehringer 00030 # </author> 00031 # 00032 # <summary> 00033 # Test suite for the CASA task specfit and tool method ia.fitprofile 00034 # </summary> 00035 # 00036 # <reviewed reviwer="" date="" tests="" demos=""> 00037 # </reviewed 00038 # 00039 # <prerequisite> 00040 # <ul> 00041 # <li> <linkto class="task_specfit.py:description">imcollapse</linkto> 00042 # </ul> 00043 # </prerequisite> 00044 # 00045 # <etymology> 00046 # Test for the specfit task and ia.fitprofile() tool method. 00047 # </etymology> 00048 # 00049 # <synopsis> 00050 # Test the specfit task and the ia.fitprofile() method upon which it is built. 00051 # </synopsis> 00052 # 00053 # <example> 00054 # 00055 # This test runs as part of the CASA python unit test suite and can be run from 00056 # the command line via eg 00057 # 00058 # `echo $CASAPATH/bin/casapy | sed -e 's$ $/$'` --nologger --log2term -c `echo $CASAPATH | awk '{print $1}'`/code/xmlcasa/scripts/regressions/admin/runUnitTest.py test_specfit[test1,test2,...] 00059 # 00060 # </example> 00061 # 00062 # <motivation> 00063 # To provide a test standard for the specfit task to ensure 00064 # coding changes do not break the associated bits 00065 # </motivation> 00066 # 00067 00068 ########################################################################### 00069 import shutil 00070 import casac 00071 from tasks import * 00072 from taskinit import * 00073 from __main__ import * 00074 import unittest 00075 import numpy 00076 00077 twogauss = "specfit_multipix_2gauss.fits" 00078 polyim = "specfit_multipix_poly_2gauss.fits" 00079 gauss_triplet = "gauss_triplet.fits" 00080 two_lorentzians = "two_lorentzians.fits" 00081 invalid_fits = "invalid_fits.im" 00082 solims = [ 00083 "amp", "ampErr", "center", "centerErr", 00084 "fwhm", "fwhmErr", "integral", "integralErr" 00085 ] 00086 birdie = "birdie.im" 00087 00088 nanvalue = 4.53345345 00089 00090 datapath=os.environ.get('CASAPATH').split()[0]+'/data/regression/unittest/specfit/' 00091 00092 def run_fitprofile ( 00093 imagename, box, region, chans, stokes, 00094 axis, mask, ngauss, poly, multifit, model="", 00095 residual="", amp="", amperr="", center="", centererr="", 00096 fwhm="", fwhmerr="", integral="", integralerr="", 00097 estimates="", logresults=True, pampest="", pcenterest="", pfwhmest="", 00098 pfix="", gmncomps=0, gmampcon="", gmcentercon="", 00099 gmfwhmcon="", gmampest=[0], gmcenterest=[0], 00100 gmfwhmest=[0], gmfix="", logfile="", pfunc="", 00101 goodamprange=[0.0], goodcenterrange=[0.0], goodfwhmrange=[0.0], 00102 sigma="", outsigma="" 00103 ): 00104 myia = iatool() 00105 myia.open(imagename) 00106 if (not myia.isopen()): 00107 myia.done() 00108 raise Exception 00109 res = myia.fitprofile( 00110 box=box, region=region, chans=chans, 00111 stokes=stokes, axis=axis, mask=mask, 00112 ngauss=ngauss, poly=poly, estimates=estimates, 00113 multifit=multifit, 00114 model=model, residual=residual, amp=amp, 00115 amperr=amperr, center=center, centererr=centererr, 00116 fwhm=fwhm, fwhmerr=fwhmerr, integral=integral, 00117 integralerr=integralerr, logresults=logresults, pampest=pampest, 00118 pcenterest=pcenterest, pfwhmest=pfwhmest, pfix=pfix, 00119 gmncomps=gmncomps, gmampcon=gmampcon, 00120 gmcentercon=gmcentercon, gmfwhmcon=gmfwhmcon, 00121 gmampest=gmampest, gmcenterest=gmcenterest, 00122 gmfwhmest=gmfwhmest, gmfix=gmfix, logfile=logfile, 00123 pfunc=pfunc, goodamprange=goodamprange, 00124 goodcenterrange=goodcenterrange, 00125 goodfwhmrange=goodfwhmrange, sigma=sigma, outsigma=outsigma ) 00126 myia.close() 00127 myia.done() 00128 return res 00129 00130 def run_specfit( 00131 imagename, box, region, chans, stokes, 00132 axis, mask, ngauss, poly, multifit, model="", 00133 residual="", amp="", amperr="", center="", centererr="", 00134 fwhm="", fwhmerr="", integral="", integralerr="", 00135 wantreturn=True, estimates="", logresults=None, 00136 pampest="", pcenterest="", pfwhmest="", pfix="", 00137 gmncomps=0, gmampcon="", gmcentercon="", 00138 gmfwhmcon="", gmampest=[0], gmcenterest=[0], 00139 gmfwhmest=[0], gmfix="", logfile="", pfunc="", 00140 goodamprange=[0.0], goodcenterrange=[0.0], goodfwhmrange=[0.0], 00141 sigma="", outsigma="" 00142 ): 00143 return specfit( 00144 imagename=imagename, box=box, region=region, 00145 chans=chans, stokes=stokes, axis=axis, mask=mask, 00146 ngauss=ngauss, poly=poly, estimates=estimates, 00147 multifit=multifit, 00148 model=model, residual=residual, amp=amp, 00149 amperr=amperr, center=center, centererr=centererr, 00150 fwhm=fwhm, fwhmerr=fwhmerr, integral=integral, 00151 integralerr=integralerr, 00152 wantreturn=wantreturn, logresults=logresults, pampest=pampest, 00153 pcenterest=pcenterest, pfwhmest=pfwhmest, pfix=pfix, 00154 gmncomps=gmncomps, gmampcon=gmampcon, 00155 gmcentercon=gmcentercon, gmfwhmcon=gmfwhmcon, 00156 gmampest=gmampest, gmcenterest=gmcenterest, 00157 gmfwhmest=gmfwhmest, gmfix=gmfix, logfile=logfile, 00158 pfunc=pfunc, 00159 goodamprange=goodamprange, 00160 goodcenterrange=goodcenterrange, 00161 goodfwhmrange=goodfwhmrange, sigma=sigma, outsigma=outsigma 00162 ) 00163 00164 class specfit_test(unittest.TestCase): 00165 00166 def setUp(self): 00167 shutil.copy(datapath + twogauss, twogauss) 00168 shutil.copy(datapath + polyim, polyim) 00169 00170 def tearDown(self): 00171 os.remove(twogauss) 00172 os.remove(polyim) 00173 00174 def checkImage(self, gotImage, expectedName): 00175 expected = iatool() 00176 expected.open(expectedName) 00177 got = iatool() 00178 if type(gotImage) == str: 00179 got.open(gotImage) 00180 else: 00181 got = gotImage 00182 self.assertTrue((got.shape() == expected.shape()).all()) 00183 gotchunk = got.getchunk() 00184 expchunk = expected.getchunk() 00185 if (numpy.isnan(gotchunk).any()): 00186 gotchunk = gotchunk.ravel() 00187 for i in range(len(gotchunk)): 00188 if isnan(gotchunk[i]): 00189 gotchunk[i] = nanvalue 00190 if (numpy.isnan(expchunk).any()): 00191 expchunk = expchunk.ravel() 00192 for i in range(len(expchunk)): 00193 if isnan(expchunk[i]): 00194 expchunk[i] = nanvalue 00195 diffData = gotchunk - expchunk 00196 self.assertTrue(abs(diffData).max() < 2e-11) 00197 self.assertTrue( 00198 ( 00199 got.getchunk(getmask=T) == expected.getchunk(getmask=T) 00200 ).all() 00201 ) 00202 gotCsys = got.coordsys() 00203 expectedCsys = expected.coordsys() 00204 diffPixels = gotCsys.referencepixel()['numeric'] - expectedCsys.referencepixel()['numeric'] 00205 self.assertTrue(abs(diffPixels).max() == 0) 00206 00207 diffRef = gotCsys.referencevalue()['numeric'] - expectedCsys.referencevalue()['numeric'] 00208 # fracDiffRef = (diffRef)/expectedCsys.referencevalue()['numeric']; 00209 self.assertTrue(abs(diffRef).max() == 0) 00210 got.close() 00211 got.done() 00212 expected.close() 00213 expected.done() 00214 00215 def test_exceptions(self): 00216 """specfit: Test various exception cases""" 00217 00218 def testit( 00219 imagename, box, region, chans, stokes, 00220 axis, mask, ngauss, poly, multifit, model, 00221 residual 00222 ): 00223 for i in [0,1]: 00224 if (i==0): 00225 self.assertRaises( 00226 Exception, run_fitprofile, imagename, 00227 box, region, chans, stokes, axis, mask, 00228 ngauss, poly, multifit, model, residual 00229 ) 00230 else: 00231 self.assertFalse( 00232 run_specfit( 00233 imagename, box, region, chans, 00234 stokes, axis, mask, ngauss, poly, 00235 multifit, model, residual 00236 ) 00237 ) 00238 # Exception if no image name given", 00239 testit( 00240 "", "", "", "", "", 2, "", False, 1, -1, "", "" 00241 ) 00242 # Exception if bogus image name given 00243 testit( 00244 "my bad", "", "", "", "", 2, "", 1, -1, False, "", "" 00245 ) 00246 # Exception if given axis is out of range 00247 testit( 00248 twogauss, "", "", "", "", 5, "", 1, -1, False, "", "" 00249 ) 00250 # Exception if bogus box string given #1 00251 testit( 00252 twogauss, "abc", "", "", "", 2, "", 1, -1, False, "", "" 00253 ) 00254 # Exception if bogus box string given #2 00255 testit( 00256 twogauss, "0,0,1000,1000", "", "", "", 2, "", 1, -1, False, "", "" 00257 ) 00258 # Exception if bogus chans string given #1 00259 testit( 00260 twogauss, "", "", "abc", "", 2, "", 1, -1, False, "", "" 00261 ) 00262 # Exception if bogus chans string given #2 00263 testit( 00264 twogauss, "", "", "0-200", "", 2, "", 1, -1, False, "", "" 00265 ) 00266 # Exception if bogus stokes string given #1 00267 testit( 00268 twogauss, "", "", "", "abc", 2, "", 1, -1, False, "", "" 00269 ) 00270 # Exception if bogus stokes string given #2 00271 testit( 00272 twogauss, "", "", "", "v", 2, "", 1, -1, False, "", "" 00273 ) 00274 # Exception if no gaussians and no polynomial specified 00275 testit( 00276 twogauss, "", "", "", "", 2, "", 0, -1, False, "", "" 00277 ) 00278 00279 def test_1(self): 00280 """Tests of averaging over a region and then fitting""" 00281 imagename = twogauss 00282 box = "" 00283 region = "" 00284 chans = "" 00285 stokes = "" 00286 axis = 2 00287 mask = "" 00288 ngauss = 2 00289 poly = -1 00290 multifit = False 00291 model = "" 00292 residual = "" 00293 for code in [run_fitprofile, run_specfit]: 00294 res = code( 00295 imagename, box, region, chans, 00296 stokes, axis, mask, ngauss, poly, 00297 multifit, model, residual 00298 ) 00299 self.assertTrue(len(res["converged"]) == 1) 00300 self.assertTrue(res["converged"][0,0,0,0]) 00301 # even though two components given, only one is fit 00302 self.assertTrue(res["ncomps"][0,0,0,0] == 1) 00303 # the fit component is a gaussian 00304 self.assertTrue(res["type"][0,0,0,0,0] == "GAUSSIAN") 00305 gs = res["gs"] 00306 self.assertAlmostEqual(gs["amp"][0,0,0,0,0], 49.7, 1, "amplitude determination failure") 00307 self.assertAlmostEqual(gs["ampErr"][0,0,0,0,0], 4.0, 1, "amplitude error determination failure") 00308 self.assertAlmostEqual(gs["center"][0,0,0,0,0], -237.7, 1, "center determination failure") 00309 self.assertAlmostEqual(gs["centerErr"][0,0,0,0,0], 1.7, 1, "center error determination failure") 00310 self.assertAlmostEqual(gs["fwhm"][0,0,0,0,0], 42.4, 1, "fwhm determination failure") 00311 self.assertAlmostEqual(gs["fwhmErr"][0,0,0,0,0], 4.0, 1, "fwhm error determination failure") 00312 00313 self.assertTrue(res["xUnit"] == "km/s") 00314 self.assertTrue(res["yUnit"] == "Jy") 00315 00316 def test_2(self): 00317 """ multipixel, two gaussian fit""" 00318 imagename = twogauss 00319 box = "" 00320 region = "" 00321 chans = "" 00322 stokes = "" 00323 axis = 2 00324 mask = "" 00325 ngauss = 2 00326 poly = -1 00327 multifit = True 00328 model = "" 00329 residual = "" 00330 for code in [run_fitprofile, run_specfit]: 00331 res = code( 00332 imagename, box, region, chans, 00333 stokes, axis, mask, ngauss, poly, 00334 multifit, model, residual 00335 ) 00336 self.assertTrue(len(res["converged"].ravel()) == 81) 00337 self.assertTrue(res["converged"].all()) 00338 self.assertTrue(res["ncomps"][0, 0, 0, 0] == 1) 00339 self.assertTrue((res["ncomps"][:, 1:, 0, 0] == 2).all()) 00340 self.assertTrue((res["ncomps"][1:, 0, 0, 0] == 2).all()) 00341 self.assertTrue((res["type"][:,:,:,:,0] == "GAUSSIAN").all()) 00342 self.assertTrue(res["type"][0, 0, 0, 0, 1] == "UNDEF") 00343 self.assertTrue((res["type"][:, 1:, 0, 0, 1] == "GAUSSIAN").all()) 00344 self.assertTrue((res["type"][1:, 0, 0, 0, 1] == "GAUSSIAN").all()) 00345 00346 self.assertTrue(res["xUnit"] == "km/s") 00347 self.assertTrue(res["yUnit"] == "Jy") 00348 00349 def test_3(self): 00350 """ Test two gaussian + one polynomial image""" 00351 imagename = polyim 00352 box = "" 00353 region = "" 00354 chans = "" 00355 stokes = "" 00356 axis = 2 00357 mask = "" 00358 ngauss = 2 00359 poly = 3 00360 multifit = True 00361 model = "" 00362 residual = "" 00363 for code in [run_fitprofile, run_specfit]: 00364 res = code( 00365 imagename, box, region, chans, 00366 stokes, axis, mask, ngauss, poly, 00367 multifit, model, residual 00368 ) 00369 self.assertTrue(len(res["converged"].ravel()) == 81) 00370 # fit #72 did not converge 00371 self.assertTrue(res["converged"][:, :7, 0, 0].all()) 00372 self.assertTrue(res["converged"][1:,8,0,0].all()) 00373 self.assertFalse(res["converged"][0, 8, 0, 0]) 00374 00375 self.assertTrue(res["xUnit"] == "km/s") 00376 self.assertTrue(res["yUnit"] == "Jy") 00377 00378 def test_4(self): 00379 """writing solution images for multipixel, two gaussian fit""" 00380 imagename = twogauss 00381 box = "" 00382 region = "" 00383 chans = "" 00384 stokes = "" 00385 axis = 2 00386 mask = "" 00387 ngauss = 2 00388 poly = -1 00389 multifit = True 00390 model = "" 00391 residual = "" 00392 [ 00393 amp, amperr, center, centererr, 00394 fwhm, fwhmerr, integral, integralerr 00395 ] = solims 00396 for code in [run_fitprofile, run_specfit]: 00397 res = code( 00398 imagename, box, region, chans, 00399 stokes, axis, mask, ngauss, poly, 00400 multifit, model, residual, amp, 00401 amperr, center, centererr, fwhm, fwhmerr, 00402 integral, integralerr 00403 ) 00404 for im in solims: 00405 self.checkImage(im, datapath + im) 00406 shutil.rmtree(im) 00407 00408 def test_5(self): 00409 """test results of multi-pixel one gaussian fit with estimates file""" 00410 imagename = twogauss 00411 ngauss=10 00412 box = "" 00413 region = "" 00414 chans = "" 00415 stokes = "" 00416 axis = 2 00417 mask = "" 00418 poly = -1 00419 estimates = datapath + "goodProfileEstimatesFormat_2.txt" 00420 multifit = True 00421 model = "" 00422 residual = "" 00423 for code in [run_fitprofile, run_specfit]: 00424 res = code( 00425 imagename, box, region, chans, 00426 stokes, axis, mask, ngauss, poly, 00427 multifit, model, residual, 00428 estimates=estimates 00429 ) 00430 # no tests yet, just confirm it runs to completion 00431 00432 def test_6(self): 00433 """test results of non-multi-pixel one gaussian fit with estimates file""" 00434 imagename = twogauss 00435 ngauss=10 00436 box = "" 00437 region = "" 00438 chans = "" 00439 stokes = "" 00440 axis = 2 00441 mask = "" 00442 poly = -1 00443 estimates = datapath + "goodProfileEstimatesFormat_2.txt" 00444 multifit = False 00445 model = "" 00446 residual = "" 00447 for code in [run_fitprofile, run_specfit]: 00448 res = code( 00449 imagename, box, region, chans, 00450 stokes, axis, mask, ngauss, poly, 00451 multifit, model, residual, 00452 estimates=estimates 00453 ) 00454 self.assertTrue(res['converged'] == 1) 00455 self.assertTrue(res['converged'][0][0][0][0]) 00456 self.assertTrue(res['ncomps'][0][0][0][0] == 1) 00457 self.assertTrue(res["type"][0,0,0,0,0] == "GAUSSIAN") 00458 gs = res["gs"] 00459 self.assertAlmostEqual(gs["amp"][0,0,0,0,0], 49.7, 1, "amplitude determination failure") 00460 self.assertAlmostEqual(gs["ampErr"][0,0,0,0,0], 4.0, 1, "amplitude error determination failure") 00461 self.assertAlmostEqual(gs["center"][0,0,0,0,0], -237.7, 1, "center determination failure") 00462 self.assertAlmostEqual(gs["centerErr"][0,0,0,0,0], 1.7, 1, "center error determination failure") 00463 self.assertAlmostEqual(gs["fwhm"][0,0,0,0,0], 42.4, 1, "fwhm determination failure") 00464 self.assertAlmostEqual(gs["fwhmErr"][0,0,0,0,0], 4.0, 1, "fwhm error determination failure") 00465 self.assertTrue(res["xUnit"] == "km/s") 00466 self.assertTrue(res["yUnit"] == "Jy") 00467 00468 def test_7(self): 00469 """test results of non-multi-pixel one gaussian fit with estimates file keeping peak fixed""" 00470 imagename = twogauss 00471 ngauss=10 00472 box = "" 00473 region = "" 00474 chans = "" 00475 stokes = "" 00476 axis = 2 00477 mask = "" 00478 poly = -1 00479 estimates = datapath + "goodProfileEstimatesFormat_3.txt" 00480 multifit = False 00481 model = "" 00482 residual = "" 00483 for code in [run_fitprofile, run_specfit]: 00484 res = code( 00485 imagename, box, region, chans, 00486 stokes, axis, mask, ngauss, poly, 00487 multifit, model, residual, 00488 estimates=estimates 00489 ) 00490 self.assertTrue(res['converged'] == 1) 00491 self.assertTrue(res['converged'][0][0][0][0]) 00492 self.assertTrue(res['ncomps'][0][0][0][0] == 1) 00493 self.assertTrue(res["type"][0,0,0,0,0] == "GAUSSIAN") 00494 self.assertAlmostEqual(res["gs"]["amp"][0,0,0,0,0], 45, 1, "amplitude determination failure") 00495 self.assertAlmostEqual(res["gs"]["ampErr"][0,0,0,0,0], 0, 1, "amplitude error determination failure") 00496 self.assertAlmostEqual(res["gs"]["center"][0,0,0,0,0], -237.7, 1, "center determination failure") 00497 self.assertAlmostEqual(res["gs"]["centerErr"][0,0,0,0,0], 1.9, 1, "center error determination failure") 00498 self.assertAlmostEqual(res["gs"]["fwhm"][0,0,0,0,0], 45.6, 1, "fwhm determination failure") 00499 self.assertAlmostEqual(res["gs"]["fwhmErr"][0,0,0,0,0], 3.8, 1, "fwhm error determination failure") 00500 self.assertTrue(res["xUnit"] == "km/s") 00501 self.assertTrue(res["yUnit"] == "Jy") 00502 00503 00504 def test_stretch(self): 00505 """specfit : test mask stretch""" 00506 imagename = twogauss 00507 yy = iatool() 00508 yy.open(imagename) 00509 mycsys = yy.coordsys().torecord() 00510 yy.done() 00511 mymask = "maskim" 00512 yy.fromshape(mymask, [9, 9, 1, 1]) 00513 yy.setcoordsys(mycsys) 00514 yy.addnoise() 00515 yy.done() 00516 yy.open(imagename) 00517 self.assertRaises( 00518 Exception, yy.fitprofile, 00519 ngauss=2, mask=mymask + ">-100", 00520 stretch=False 00521 ) 00522 zz = specfit( 00523 imagename, ngauss=2, mask=mymask + ">-100", 00524 stretch=False 00525 ) 00526 self.assertTrue(zz == None) 00527 zz = yy.fitprofile( 00528 ngauss=2, mask=mymask + ">-100", 00529 stretch=True 00530 ) 00531 self.assertTrue(len(zz.keys()) > 0) 00532 yy.done() 00533 zz = specfit( 00534 imagename, ngauss=2, mask=mymask + ">-100", 00535 stretch=True 00536 ) 00537 self.assertTrue(len(zz.keys()) > 0) 00538 00539 def test_8(self): 00540 """ Test two gaussian + one polynomial image with estimates""" 00541 00542 estimates = datapath + "poly+2gauss_estimates.txt" 00543 00544 for code in [run_fitprofile, run_specfit]: 00545 res = code( 00546 imagename=polyim, box="", region="", chans="", 00547 stokes="", axis=2, mask="", ngauss=0, poly=3, 00548 multifit=True, model="", residual="", estimates=estimates 00549 ) 00550 self.assertTrue(len(res["converged"].ravel()) == 81) 00551 # fit #0 did not converge 00552 self.assertTrue(res["converged"][1:, :, 0, 0].all()) 00553 self.assertTrue(res["converged"][0,1:,0,0].all()) 00554 self.assertFalse(res["converged"][0, 0, 0, 0]) 00555 00556 self.assertTrue(res["xUnit"] == "km/s") 00557 self.assertTrue(res["yUnit"] == "Jy") 00558 gs = res["gs"] 00559 center = gs["center"] 00560 center[0,0,0,0,0] = nanvalue 00561 center[0,0,0,0,1] = nanvalue 00562 00563 amp = gs["amp"] 00564 amp[0,0,0,0,0] = nanvalue 00565 amp[0,0,0,0,1] = nanvalue 00566 00567 fwhm = gs["fwhm"] 00568 fwhm[0,0,0,0,0] = nanvalue 00569 fwhm[0,0,0,0,1] = nanvalue 00570 00571 for code in [run_fitprofile, run_specfit]: 00572 res = code( 00573 imagename=polyim, box="", region="", chans="", 00574 stokes="", axis=2, mask="", ngauss=0, poly=3, 00575 multifit=True, model="", residual="", estimates="", 00576 pampest=[50, 10], pcenterest=[90, 30], pfwhmest=[10, 7] 00577 ) 00578 self.assertTrue(len(res["converged"].ravel()) == 81) 00579 # fit #0 did not converge 00580 self.assertTrue(res["converged"][1:, :, 0, 0].all()) 00581 self.assertTrue(res["converged"][0,1:,0,0].all()) 00582 self.assertFalse(res["converged"][0, 0, 0, 0]) 00583 00584 self.assertTrue(res["xUnit"] == "km/s") 00585 self.assertTrue(res["yUnit"] == "Jy") 00586 gs = res["gs"] 00587 gs["center"][0,0,0,0,0] = nanvalue 00588 gs["center"][0,0,0,0,1] = nanvalue 00589 self.assertTrue((gs["center"] == center).all()) 00590 00591 gs["amp"][0,0,0,0,0] = nanvalue 00592 gs["amp"][0,0,0,0,1] = nanvalue 00593 self.assertTrue((abs(gs["amp"] - amp) < 1e-13).all()) 00594 00595 gs["fwhm"][0,0,0,0,0] = nanvalue 00596 gs["fwhm"][0,0,0,0,1] = nanvalue 00597 self.assertTrue((abs(gs["fwhm"] - fwhm) < 1e-13).all()) 00598 00599 def test_9(self): 00600 """Polynomial fitting, moved from imagetest_regression.py""" 00601 shape = [16,16,128] 00602 imname = 'ia.fromshape.image' 00603 myim = ia.newimagefromshape(shape=shape) 00604 myim.set(pixels='1.0') # 00605 residname = 'ia.fromshape.resid' 00606 fitname = 'ia.fromshape.fit' 00607 res = myim.fitprofile (multifit=True, residual=residname, model=fitname, poly=0, ngauss=0, axis=2) 00608 myim.done() 00609 myim.open(residname) 00610 pixels = myim.getchunk() 00611 self.assertTrue(pixels.shape == tuple(shape)) 00612 self.assertTrue((pixels == 0).all()) 00613 myim.done() 00614 00615 def test_10(self): 00616 """test results of non-multi-fit gaussian triplet""" 00617 imagename=datapath+gauss_triplet 00618 gmampcon = [0.7, 0.55] 00619 gmcentercon = [52, 0] 00620 for code in [run_fitprofile, run_specfit]: 00621 res = code( 00622 imagename=imagename, box="", region="", chans="", 00623 stokes="", axis=2, mask="", ngauss=0, poly=-1, 00624 multifit=False, model="", residual="", estimates="", 00625 gmncomps=3, gmampest=[1.2, 0.8, 0.6], 00626 gmcenterest=[20, 72, 100], gmfwhmest=[4, 4, 4], 00627 gmampcon=gmampcon, gmcentercon=gmcentercon 00628 ) 00629 self.assertTrue(res["type"].ravel() == "GAUSSIAN MULTIPLET") 00630 gm0 = res["gm0"] 00631 exp = [4.15849, 2.91095, 2.28717] 00632 got = gm0["amp"].ravel() 00633 for i in [0, 1, 2]: 00634 self.assertAlmostEqual(got[i], exp[i], 5) 00635 exp = [1149.73, 1138.76, 1133.66] 00636 got = gm0["center"].ravel() 00637 for i in [0, 1, 2]: 00638 self.assertAlmostEqual(got[i], exp[i], 2) 00639 exp = [5.75308, 4.09405, 3.93497] 00640 got = gm0["fwhm"].ravel() 00641 for i in [0, 1, 2]: 00642 self.assertAlmostEqual(got[i], exp[i], 5) 00643 exp = [0.0301945, 0.0211362, 0.016607] 00644 got = gm0["ampErr"].ravel() 00645 for i in [0, 1, 2]: 00646 self.assertAlmostEqual(got[i], exp[i], 7) 00647 exp = [0.0221435, 0.0221435, 0.0475916] 00648 got = gm0["centerErr"].ravel() 00649 for i in [0, 1, 2]: 00650 self.assertAlmostEqual(got[i], exp[i], 7) 00651 exp = [0.0556099, 0.085414, 0.0987483] 00652 got = gm0["fwhmErr"].ravel() 00653 for i in [0, 1, 2]: 00654 self.assertAlmostEqual(got[i], exp[i], 7) 00655 self.assertAlmostEqual( 00656 gm0["amp"].ravel()[0]*gmampcon[0], 00657 gm0["amp"].ravel()[1], 7 00658 ) 00659 self.assertAlmostEqual( 00660 gm0["ampErr"].ravel()[0]*gmampcon[0], 00661 gm0["ampErr"].ravel()[1], 7 00662 ) 00663 self.assertAlmostEqual( 00664 gm0["amp"].ravel()[0]*gmampcon[1], 00665 gm0["amp"].ravel()[2], 7 00666 ) 00667 self.assertAlmostEqual( 00668 gm0["ampErr"].ravel()[0]*gmampcon[1], 00669 gm0["ampErr"].ravel()[2], 7 00670 ) 00671 ia.open(imagename) 00672 mc = ia.coordsys() 00673 ia.done() 00674 restfreq = mc.restfrequency()["value"][0] 00675 dv = mc.increment()["numeric"][2] 00676 increment = -dv/restfreq*299797 00677 self.assertAlmostEqual( 00678 gm0["center"].ravel()[0] + gmcentercon[0]*increment, 00679 gm0["center"].ravel()[1], 3 00680 ) 00681 self.assertAlmostEqual( 00682 gm0["centerErr"].ravel()[0], 00683 gm0["centerErr"].ravel()[1], 7 00684 ) 00685 00686 def test_11(self): 00687 """test results of multi-fit gaussian triplet""" 00688 imagename=datapath+gauss_triplet 00689 gmampcon = [0.7, 0.55] 00690 gmcentercon = [52, 0] 00691 logfile = "mylog.txt" 00692 i = 1 00693 for code in [run_fitprofile, run_specfit]: 00694 res = code( 00695 imagename=imagename, box="", region="", chans="", 00696 stokes="", axis=2, mask="", ngauss=0, poly=-1, 00697 multifit=True, center="center", 00698 centererr="centerErr", fwhm="fwhm", 00699 fwhmerr="fwhmErr", amp="amp", amperr="ampErr", 00700 integral="integral", integralerr="integralErr", 00701 gmncomps=3, gmampest=[1.2, 0.1, 0.1], 00702 gmcenterest=[20, 0, 100], gmfwhmest=[4, 4, 4], 00703 gmampcon=gmampcon, gmcentercon=gmcentercon, 00704 logfile=logfile 00705 ) 00706 for image in ( 00707 "center", "centerErr", "fwhm", "fwhmErr", "amp", 00708 "ampErr", "integral", "integralErr" 00709 ): 00710 self.checkImage( 00711 image + "_gm", datapath + image + "_gm" 00712 ) 00713 # appending, second time through size should double 00714 self.assertTrue(os.path.getsize(logfile) > 3e4*i) 00715 i = i+1 00716 00717 def test_12(self): 00718 """test results of lorentzian fitting""" 00719 imagename=datapath+two_lorentzians 00720 pamp = [1, 7] 00721 pcen = [30, 111] 00722 pfwhm = [4, 4] 00723 pfunc = ["l", "l"] 00724 00725 logfile = "two_lorentzian_fit.log" 00726 i = 1 00727 for code in [run_fitprofile, run_specfit]: 00728 res = code( 00729 imagename=imagename, box="", region="", chans="", 00730 stokes="", axis=2, mask="", ngauss=0, poly=-1, 00731 multifit=True, center="center", 00732 centererr="centerErr", fwhm="fwhm", 00733 fwhmerr="fwhmErr", amp="amp", amperr="ampErr", 00734 integral="integral", integralerr="integralErr", 00735 pampest=pamp, pcenterest=pcen, pfwhmest=pfwhm, 00736 logfile=logfile, pfunc=pfunc 00737 ) 00738 for image in ( 00739 "center", "centerErr", "fwhm", "fwhmErr", "amp", 00740 "ampErr", "integral", "integralErr" 00741 ): 00742 self.checkImage( 00743 image + "_ls", datapath + image + "_ls" 00744 ) 00745 # appending, second time through size should double 00746 self.assertTrue(os.path.getsize(logfile) > 2e4*i) 00747 i = i+1 00748 00749 def test_13(self): 00750 """test setting solution parameter validities """ 00751 imagename = datapath + invalid_fits 00752 i = 0 00753 goodfwhmrange= [[0], [2,20]] 00754 for gfr in goodfwhmrange: 00755 for code in [run_fitprofile, run_specfit]: 00756 res = code( 00757 imagename=imagename, box="", region="", chans="", 00758 stokes="", axis=3, mask="", ngauss=2, poly=-1, 00759 multifit=True, logresults=False, goodfwhmrange=gfr 00760 ) 00761 if i == 0: 00762 exp = 16 00763 else: 00764 exp = 6 00765 #print "valid " + str(res["valid"]) 00766 #print "*** fwhm " + str(res["gs"]["fwhm"]) 00767 self.assertTrue(res["valid"].sum() == exp) 00768 i = i+1 00769 00770 def test_14(self): 00771 imagename = datapath + birdie 00772 sigmaimage = "sigma.im" 00773 myia = iatool() 00774 myia.open(imagename) 00775 bb = myia.getchunk() 00776 myia.done() 00777 fullsigma = bb 00778 fullsigma[:] = 1 00779 s = [] 00780 s.append(fullsigma[:]) 00781 s.append(fullsigma[0, 0, :, 0].reshape([1, 1, 100, 1])) 00782 s.append(fullsigma[0, 0, :, 0]) 00783 outsigma = "outsigma.im" 00784 expfwhmtrue = { 00785 "0": [ 61.77550675, 41.74020775], 00786 "0.1" : [82.55515563, 36.4484157], 00787 "1": [ 61.89796905, 41.69283371], 00788 "100": [ 61.77551883, 41.74020306] 00789 } 00790 expfwhmerrtrue = { 00791 "0": [ 2.88886469e-07, 4.18990637e-07], 00792 "0.1": [ 1.96169836, 2.99212108], 00793 "1": [ 0.30788521, 0.44558165], 00794 "100": [ 0.00307781, 0.00446394] 00795 } 00796 expfwhmfalse = { 00797 "0": 42.43941294, 00798 "0.1": 42.43941294, 00799 "1": 42.43941294, 00800 "100": 42.43941294 00801 } 00802 expfwhmerrfalse = { 00803 "0": 4.0707913, 00804 "0.1": 6.09143172, 00805 "1": 4.07523664, 00806 "100": 4.04975604 00807 } 00808 for sigma in (s): 00809 for birdiesigma in [0, 0.1, 1, 100]: 00810 fullsigma[:, :, 50, :] = birdiesigma 00811 mymax = max(1, birdiesigma) 00812 for i in range(2): 00813 sig = sigma 00814 if (sig.ndim == 1): 00815 sig[50] = birdiesigma 00816 else: 00817 sig[:,:,50,:] = birdiesigma 00818 if i == 1: 00819 myia.fromarray(sigmaimage, sig, overwrite=T) 00820 myia.done() 00821 sig = sigmaimage 00822 for code in [run_fitprofile, run_specfit]: 00823 for multifit in [True, False]: 00824 res = code( 00825 imagename=imagename, box="", region="", chans="", 00826 stokes="", axis=2, mask=imagename + "<1000", ngauss=2, poly=-1, 00827 multifit=multifit, logresults=False, sigma=sig, 00828 outsigma=outsigma 00829 ) 00830 if multifit: 00831 self.assertTrue((abs(res["gs"]["fwhm"][8,8,0,0,:] - expfwhmtrue[str(birdiesigma)]) < 1e-7).all()) 00832 self.assertTrue((abs(res["gs"]["fwhmErr"][8,8,0,0,:] - expfwhmerrtrue[str(birdiesigma)]) < 1e-7).all()) 00833 else: 00834 self.assertTrue(abs(res["gs"]["fwhm"][0,0,0,0,0] - expfwhmfalse[str(birdiesigma)]) < 1e-7) 00835 self.assertTrue(abs(res["gs"]["fwhmErr"][0,0,0,0,0] - expfwhmerrfalse[str(birdiesigma)]) < 1e-7) 00836 myia.open(outsigma) 00837 if (birdiesigma == 0 or birdiesigma == 1): 00838 self.assertTrue((mymax*myia.getchunk() == fullsigma).all()) 00839 else: 00840 self.assertTrue(((mymax*myia.getchunk() - fullsigma)/fullsigma < 1e-7).all()) 00841 myia.remove() 00842 00843 00844 00845 00846 def suite(): 00847 return [specfit_test]
1.8.0
